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Abstract 


Convergent disk migration has long been suspected to be responsible for forming planetary systems with a chain of 
mean-motion resonances (MMRs). Dynamical evolution over time could disrupt the delicate resonant 
configuration. We present TOI-1136, a 700 + 150 Myr old С star hosting at least six transiting planets between 
~2 and 5 Ro. The orbital period ratios deviate from exact commensurability by only 1074, smaller than the ~1072 
deviations seen in typical Kepler near-resonant systems. A transit-timing analysis measured the masses of the 
planets (3-8М-) and demonstrated that the planets in TOI-1136 are in true resonances with librating resonant 
angles. Based on a Rossiter-McLaughlin measurement of planet d, the star's rotation appears to be aligned with the 
planetary orbital planes. The well-aligned planetary system and the lack of a detected binary companion together 
suggest that TOI-1136's resonant chain formed in an isolated, quiescent disk with no stellar flyby, disk warp, or 
significant axial asymmetry. With period ratios near 3:2, 2:1, 3:2, 7:5, and 3:2, TOI-1136 is the first known 
resonant chain involving a second-order MMR (7:5) between two first-order MMRs. The formation of the delicate 
7:5 resonance places strong constraints on the system's migration history. Short-scale (starting from ~0.1 au) 
Type-I migration with an inner disk edge is most consistent with the formation of TOI-1136. A low disk surface 
density (%1 au Í 107g cm 2; lower than the minimum-mass solar nebula) and the resultant slower migration rate 
likely facilitated the formation of the 7:5 second-order MMR. 
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Unified Astronomy Thesaurus concepts: Exoplanet dynamics (490); Exoplanet evolution (491); Exoplanet 


migration (2205); Exoplanet formation (492) 


1. Introduction 


Disk migration is predicted to be a common stage of planet 
formation: in most scenarios the net effect is migration toward 
the central star (Goldreich & Tremaine 1979; Lin & 
Papaloizou 1986; Ward 1997; McNeil et al. 2005; Terquem 
& Papaloizou 2007; Nelson 2018). A pair of planets may 
become locked into a mean-motion resonance (MMR) if the 
migration is slow (adiabatic) and convergent (outer planets 
catching up with the inner planet) This process can be 
extended to capture multiple planets in a chain of resonance 
(see Kley2012, and references therein). Different studies using 
adiabatic perturbation theory (Henrard 1982; Batygin 2015), 
modified N-body integration (e.g., Lee & Peale 2002; Terquem 
& Papaloizou 2007), and hydrodynamic simulations (e.g., Kley 
et al. 2005; McNeil et al. 2005; Cresswell & Nelson 2008; 
Ogihara & Ida 2009; Ataiee & Kley 2020) all came to the same 
conclusion that convergent disk migration consistently gen- 
erates compact, first-order resonant chains of planets. This 
process of resonant capture is considered to be so effective and 
robust that it is difficult to understand why only a few percent 
of Kepler multiplanet systems are near first-order MMR 
(Fabrycky et al. 2014). Upon closer examination, most of 
these systems still show 196—296 positive deviation from 
perfect period commensurability. Transit-timing-variation 
(TTV) modeling (e.g., Hadden & Lithwick 2017) has shown 
that most of these systems are near-resonant (with circulating 
resonant angles) rather than being truly resonant (librating 
resonant angles). 

Planetesimal scattering (Chatterjee & Ford 2015), tidal 
dissipation (Lithwick & Wu 2012; Batygin & Morbidelli 
2013a), secular chaos (Petrovich et al. 2018), and orbital 
instability (Pu & Wu 2015; Izidoro et al. 2017; Goldberg & 
Batygin 2022) are some of the possible mechanisms for breaking 
migration-induced resonances as planetary systems mature. Some 
of these processes may take as long as billions of years to 
manifest. One might expect, therefore, that when the Kepler 
multiplanet systems were younger, they were also closer to 
resonance or truly resonant. In this paper, we present a young 
system that is deep in resonance (the observed orbital period ratios 
are close to small integer ratios; the relevant resonant angles are 
also librating). TOI-1136 has a resonant chain of at least six 


transiting planets, all of which display TTVs. The planets’ orbital 
period ratios deviate from the perfect integer period ratio by 1074. 
With an age of only 700 Myr, TOI-1136 may still record a pristine 
orbital architecture produced by convergent disk migration, before 
subsequent dynamical evolution has had the chance to disrupt the 
resonance. We present in this paper a series of observations and 
dynamical modeling to characterize the system and explore how 
the system formed and dynamically evolved. 

The paper is organized as follows. Section 2 characterizes 
the host star, establishes its youth, and puts limits of the 
presence of a stellar companion. Section 3 presents a Rossiter— 
McLaughlin measurement of planet d. Sections 4 and 5 contain 
our analyses of the transit signal and transit-timing variations. 
Section 6 describes a series of dynamical models to investigate 
the dynamical stability, resonant configuration, disk migration, 
and resonant repulsion of TOI-1136. Section 7 discusses the 
implications for the formation and evolution of TOI-1136 in 
relation to other multiplanet systems. Section 8 is a brief 
summary of the paper. 


2. Host Star Properties 
2.1. Spectroscopic and Stellar Parameters 


We obtained three high-resolution, high-signal-to-noise-ratio 
(S/N), iodine-free spectra of TOI-1136 with the High 
Resolution Echelle Spectrometer (HIRES) on the 10m Keck 
I telescope (Vogt et al. 1994). We employed SpecMatch- 
Зуп (for details see Petigura et al. 2017) to extract the 
spectroscopic parameters (Тең, logg, and [Fe/H]) of the host 
star. The results are listed in Table 1. The cross-correlation 
function of our HIRES spectra ruled out a spectroscopic binary 
that contributes more than 1% of the observed flux. 

To derive the stellar parameters, including the mass and 
radius of the host star, we fitted the measured spectroscopic 
parameters with Gaia parallax information (Gaia Collaboration 
et al. 2018) in the Lsoclassify package (Huber et al. 2017). 
Our procedure was similar to that presented in Fulton & 
Petigura (2018). We summarize the stellar parameters in 
Table 1. Tayar et al. (2022) showed that between different 
theoretical model grids, the systematic uncertainties from 


23 https: //github.com/petigura/specmatch-syn 
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Table 1 
Stellar Parameters of TOI-1136 


Value and 68.3% Cred- 


Parameters ible Interval Reference 
TIC ID 142276270 A 
R.A. 12:48:44.38 A 
Decl. +64:51:18.99 A 
V (mag) 9.534 + 0.003 A 
K (mag) 8.034 + 0.021 A 
Distance (pc) 84.5362 + 0.158 A 
Effective Temperature Т, (К) 5770 + 50 B 
Surface Gravity log g (dex) 4.47 + 0.04 B 
Iron Abundance [Fe/H] (dex) 0.07 + 0.06 B 
Rotational Broadening v sin i, 6.7 + 0.6 B 
(km s^!) 
Stellar Radius А, ( Ro) 0.968 + 0.036 B 
Stellar Mass M, (Mo) 1.022 + 0.027 B 
Stellar Density р, (р) 1.11 + 0.12 B 
Limb Darkening 41 0.38 + 0.16 B 
Limb Darkening 4> 0.24 + 0.11 B 
Activity Indicator Suy 0.32 + 0.03 B 
Activity Indicator logRfik -4.49 + 0.05 B 
Age (Myr) from Gyrochronology, 700 + 100 B 


Activity Indicator, and Lithium 


Note. A:TICv8 (Stassun et al. 2019); B: this work. 


Isoclassify could potentially amount to ~2% in Тең, ~4% 
in M,, and ~5% in R,. We caution the readers that these 
systematic uncertainties are not explicitly included in Table 1. 


2.2. Rotation Period 


We measured the rotation period of TOI-1136 from the 
rotational modulation seen in the Transiting Exoplanet Survey 
Satellite (TESS) light curve. With a Lomb-Scargle period- 
ogram (Lomb 1976; Scargle 1982), we measured a period of 
Pro = 8.7 + 0.1 days for the strongest peak in the periodogram. 
The corresponding flux variation has an amplitude of about 1% 
(see Figure 22). 

We estimated the age of the system using gyrochronology. 
Given a 8.7 + 0.1 day rotation period for a star like TOI-1136, 
the gyrochronal relation from Schlaufman (2010) yields an age 
of 610 + 15 Myr. Alternatively, if one follows Mamajek & 
Hillenbrand (2008), the estimated age is 700 + 20 Myr. То 
leverage the latest empirical results, we put TOI-1136 on a 
rotation versus de-reddened color diagram to compare against 
young clusters with precise rotation period measurements 
(Figure 1). Given Свр — Ggp = 0.81 and ignoring reddening 
due to the ^85 pc distance, TOI-1136 rotates at roughly the 
same rate as stars with comparable color in Praesepe (670 Myr 
old; Douglas et al. 2017). It rotates slower than any comparable 
star in M48 (450 Myr old; Barnes et al. 2015), and faster than 
any comparable star in NGC 6811 (1 Gyr old; Curtis et al. 
2019). Given TOI-1136's overlap with the stars in the Praesepe 
cluster, we conclude that the age of TOI-1136 is 2700 Myr. 
This estimate is tied to Praesepe's age, which could be as high 
as 800 Myr (Brandt & Huang 2015). 


2.3. Lithium Absorption 


We modeled the lithium absorption in our Keck/HIRES 
spectra of TOI-1136 to corroborate the youth of system. We 
modeled the Li I doublet at 6708 А as well as the nearby Fe I 
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Figure 1. Rotation period and de-reddened Gaia Ggp — Сер color of TOI-1136 
(yellow star) and stars within selected young clusters. Based on this 
"gyrochronal" comparison, TOI-1136 is likely younger than NGC 6811 
(1 Gyr old; Curtis et al. 2019) and older than M48 (450 Myr old; Barnes 
et al. 2015). The current rotation period of TOI-1136 (8.7 + 0.1 day) suggests 
an age similar or slightly older than Praesepe (670 Myr old; Douglas 
et al. 2017). 


line simultaneously. Following the procedure of Bouma et al. 
(2021), we estimated an equivalent width (EW) of 67.9 + 1.0 
mA. This Li EW is again consistent with the Praesepe cluster 
(see Figure 7 of Bouma et al. 2021) and is higher than that of 
most field stars (see also Figure 5 of Berger et al. 2018). 


2.4. Ca HK Emission and Adopted Age 


Chromospheric emission lines can provide further constraint 
on the youth of TOI-1136. We analyzed the Ca II H&K lines in 
our HIRES spectra and extracted the Sy and log Күк values 
using the method of Isaacson & Fischer (2010). TOI-1136 has 
enhanced stellar activity compared to field stars: we obtained a 
mean Syk = 0.32 + 0.03 and log Rik = —4.49 + 0.05 (field 
stars of similar spectral type typically have log Rix = — 5.0; 
Isaacson & Fischer 2010). We converted the log Кук to an 
estimate of the age of the host star. We followed the relation 
linking B — V color, log Аук, and age calibrated by Mamajek & 
Hillenbrand (2008). The age of TOI-1136 was estimated to be 
570 + 200 Myr, consistent with the age from gyrochronology 
and Li absorption. We combined the various age indicators by 
taking a weighted average, and we enlarged the formal 
uncertainty to reflect the systematic uncertainties in the 
different methods to arrive at an age for TOI-1136 of 
700 + 150 Myr. 


2.5. Cluster Membership 


Given its youth, TOI-1136 may be part of a young comoving 
group. We checked the proper motion of TOI-1136 for 
comoving groups against Banyan-» (Gagne et al. 2018) as 
well as the more recent compilation of open clusters and 
moving groups by Bouma et al. (2022). No match was found. 
We also used the Python package COMOVE (Tofflemire et al. 
2021) to search for comoving stars. We limited our search to a 
radius of 25 pc in spatial separation. COMOVE returned 11 stars 
with a tangential velocity difference <2 kms ! within this 
search volume. The closest had a 3D separation of about 17 pc. 
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Figure 2. Contrast curve as function of the radial separation for TOI-1136 
using the AO imaging from Gemini/NIRI in the K band. No stellar companion 
was identified. 


These separations could not establish a firm kinematic 
connection between these stars and TOI-1136. 


2.6. High-resolution Imaging 


To rule out nearby stellar companion, we performed a series of 
high-resolution imaging on TOI-1136 (see the Appendix). We 
highlight here the Adaptive Optics (AO) imaging observation on 
Gemini Near-Infrared Imager (NIRI; Hodapp et al. 2003) on UT 
2019 December 6. We obtained nine frames, each with an 
exposure time of 1.8 s in the Bry band. We dithered the frames by 
2" in a 2D grid. The data were reduced with a custom IDL routine 
that removes bad pixels, subtracts the sky background, flattens the 
field, and co-adds the frames. No stellar companion was seen 
anywhere in the combined image (total field of view of 
7-26" x 26"). We also performed an injection/recovery test to 
quantify the sensitivity of the AO observation. The resultant 
sensitivity curve is shown in Figure 2. We can rule out 
companions with Amag of 6.4 at separations larger than 0"5. 


2.7. A Single Star 


TOI-1136 has no reported a visual or comoving companion 
on SIMBAD, VIZIER, or Gaia data release 3 (DR3; Kervella 
et al. 2022). Gaia DR3 astrometry provides additional 
information on the possibility of inner companions that may 
have gone undetected by either Gaia or the high-resolution 
imaging. The Gaia renormalized unit weight error (RUWE) is a 
metric, similar to a reduced chi-square, where values that are 
<1.4 indicate that the Gaia astrometric solution is consistent 
with the star being single, whereas КОМЕ values 21.4 may 
indicate an astrometric excess noise, possibly caused the 
presence of an unseen companion (e.g., Ziegler et al. 2020). 
TOI-1136 has a Gaia DR3 RUWE value of 0.99 indicating that 
the astrometric fit is consistent with a single-star model. The 
lack of a spectroscopic (spectra), blended (AO), visual 
(SIMBAD), and comoving (Gaia) companion indicate that 
TOI-1136 is likely a single star. 


3. Rossiter-McLaughlin Observation 


TOL1136 is amenable to a Rossiter-McLaughlin (RM) 
measurement given its large rotational broadening, bright 
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V-band magnitude, and relatively long transit duration. More- 
over, it provides a rare chance to obtain a stellar obliquity 
measurement for a young planetary system with a resonant 
chain of planets. We observed a total of 52 spectra of TOI-1136 
with the HIRES on the 10 m Keck I telescope (Keck/HIRES; 
Vogt et al. 1994) on the night of UTC 2022 March 11 during a 
transit of TOI-1136 d as part of the TESS Keck Survey (TKS; 
see Chontos et al. 2022). We obtained the spectra with the 
iodine cell in the light path. The dense and well-measured 
molecular lines serve to anchor the wavelength solution and the 
model of the line spread function. Each exposure lasted about 
500 s and reached a median S/N of 200 per reduced pixel near 
5500 À. We had previously obtained a series of iodine-free 
spectra, which were used to create a high-S/N template stellar 
spectrum for radial velocity extraction. The radial velocities 
were extracted using our standard HIRES forward-modeling 
pipeline (Howard et al. 2010). The extracted RVs and 
uncertainties are reported in Table 5. 

We used our best-fit transit model from the TESS light 
curves (Section 4) to assist in the modeling of the RM effect. 
Specifically, we modeled the phase-folded, transit-timing- 
variation (TTV)-adjusted TESS transits of planet d simulta- 
neously with the RM effect. The model for the RM effect 
included the time of conjunction as a free parameter to account 
for the large TTVs. Reassuringly, the best-fit mid-transit time 
of the RM measurement confirmed the TTV of planet d and 
followed the trend that we expected from the TESS data. Our 
RM model closely follows the prescription of Hirano et al. 
(2011). In addition to the usual transit parameters modeled in 
Section 4, the RM model also requires the following 
parameters: the sky-projected obliquity A, the projected 
rotational velocity v sini,, a linear function of time to describe 
the local radial velocity (RV) trend with an offset y, and the 
local gradient у. An RV jitter term was also included to 
subsume any additional astrophysical or instrumental noise. No 
clear sign of a red noise component was seen in the RM 
residuals (Figure 3); we therefore adopted а simple ж 
likelihood function with a penalty term for the jitter parameter 
(e.g., Howard et al. 2013). We found the best-fit model using 
the Levenberg-Marquardt method implemented in the 
Python package 1mfit (Newville et al. 2014). 

To sample the posterior distribution, we used the Markov 
Chain Monte Carlo (MCMC) technique implemented in emcee 
(Foreman-Mackey et al. 2013). We launched 128 walkers near 
the best-fit model, and ran them for 10,000 links. We used the 
Gelman-Rubin convergence statistic (Gelman et al. 2014) to 
assess convergence of our MCMC process. The statistic was 
below 1.01 for each parameter by the end of the process, 
indicating good convergence. The results are summarized in 
Table 10. In short, we found that TOI-1136 d has a sky-projected 
obliquity А of 5? + 5°, consistent with zero. Moreover, the RM 
modeling provided a consistent but tighter constraint on the 
rotational broadening v sin i, —6.74-0.6 km s7" compared to the 
spectroscopic value (у sini, = 5.3 + 1.3 kms '). Combining 
vsini,, the stellar radius, and the stellar rotation period from 
TESS, we placed a constraint on the stellar inclination sini, 
(Masuda & Winn 2020). Following the procedure outlined by 
Albrecht et al. (2021), we found that the stellar obliquity Y is 
consistent with being zero, with an upper limit of 28? at a 9596 
credible level. We also performed an independent RM 
measurement of TOI-1136 d on High Accuracy Radial velocity 
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Figure 3. The measured Rossiter-McLaughlin effect of TOI-1136 d suggests a well-aligned orbit with a sky-projected obliquity of A = 5° + 5°. Taking the stellar 
rotation period, the stellar radius, and v sini, measurements into account, the stellar obliquity of TOI-1136 d is consistent with being 0° with a 95% upper limit of 28°. 
The black points are our HIRES measurements. The red curve is the best-fit model; blue curves are random posterior draws. The mid-transit times from the RM 
measurement confirmed and followed the trend of the TTV seen in the TESS data (Figure 5). 


Planetary Searcher North (HARPS-N) that yielded a consistent 
result. The details are outlined in the Appendix. 


4. TESS Observations 


TOI-1136 (TIC 142276270) was observed by the TESS 
mission (Ricker et al. 2014) in Sectors 14, 15, 21, 22, 41, and 
48 from UT Jul 18 2019 to Feb 25 2022. Our analysis was 
based on the 2 min cadence light curve reduced by the TESS 
Science Processing Operations Center (SPOC; Jenkins et al. 
2016) available on the Mikulski Archive for Space Telescopes 
website. It can be accessed on DOI: 10.17909 /t9-nmc8-f686. 
We experimented with both the Simple Aperture Photometry 
(SAP; Twicken et al. 2010; Morris et al. 2020) and the 
Presearch Data Conditioning Simple Aperture Photometry 
(PDCSAP; Smith et al. 2012; Stumpe et al. 2012, 2014) 
versions of the light curves. We chose to present the results 
based on the SAP light curve in this paper. The SAP light curve 
preserves the stellar variability much better, while both 
versions produced nearly identical transit fits. We minimized 
the influence of anomalous data by excluding cadences with 
nonzero Quality flags. 


4.1. Transit Modeling 


The TESS team reported four transiting planet candidates 
(Guerrero et al. 2021) with orbital periods of 6.3 (TOI- 
1136.02), 12.5 (TOI-1136.01), 18.8 (TOI-1136.04), and 26.3 


^? https: / [ archive.stsci.edu 


(TOI-1136.03) days, based on the Threshold Crossing Events 
produced in the SPOC transit search (e.g., Jenkins et al. 2020). 
The ExoFOP website^? reported two additional planets on 4.2 
and 39.5 day orbits identified by the community (ExoFOP 
website). We confirmed the detection of these candidates with 
an independent box-least-square search (BLS; Kovacs et al. 
2002) previously used in Dai et al. (2021). 

We realized that TOI-1136 may display large TTVs given 
how close the orbital periods are to resonance (see Section 7.1). 
We employed the Python package Batman (Kreidberg 2015) 
to model the transit light curves. The precise stellar density 
derived in Section 2 served as a prior in our transit modeling. A 
precise stellar density prior assists transit modeling by 
mitigating the degeneracy in the semimajor axis, impact 
parameter, and orbital eccentricity (Seager & Mallen-Ornelas 
2003). We adopted a quadratic limb-darkening profile in the 
reparameterization of д; and 42 by Kipping (2013) for efficient 
sampling. We imposed a Gaussian prior (width — 0.3) on the 
limb-darkening coefficients centered on the theoretical values 
from EXOFAST (Eastman et al. 2013). The mean stellar density 
and the limb-darkening coefficients are the three global 
parameters shared by all planets in TOI-1136. Each planet 
has its usual transits parameters: the orbital period Porp, the 
time of conjunction Т, the planet-to-star radius ratio Rp /R,, the 
scaled orbital distance a/R,, the transit impact parameter b, the 
orbital eccentricity e, and the argument of the pericenter w. 


20 https: / /exofop.ipac.caltech.edu 
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The first step in our transit modeling was to remove any 
stellar variability and instrumental flux variation by fitting a 
cubic spline of length 0.5 day to the TESS light curve. Before 
fitting the spline, we removed any data points within 2 times 
the transit duration T,4 around each transit (and TTVs were 
accounted for in subsequent iterations of this process). The 
original light curve (with transits) was then divided by the 
spline fit. Figure 22 shows the original TESS light curve, the 
spline fit, and the detrended light curve. Visual inspection 
confirmed that the detrending procedure was successful, with 
no obvious distortions of the transit light curve. 

The next step was to fit the transits of each planet assuming a 
constant orbital period. We obtained the best-fit model with the 
Levenberg-Marquardt method implemented іп the 
Python package 1mfit (Newville et al. 2014). The best-fit 
model served as a template when we fitted for the mid-transit 
time of each individual transit. During the fit for each transit, 
the only free parameters were the mid-transit time and three 
parameters of a quadratic function of time that accounts for any 
residual out-of-transit flux variations. In TOI-1136, there are 
often cases where transits of different planets partially overlap 
with each other. In those cases, we fitted the involved planets 
simultaneously. The loss of light was assumed to be the sum of 
the losses due to each planet, without accounting for possible 
planet-planet eclipses (e.g., Hirano et al. 2012). We delay a 
thorough investigation of possible planet-planet eclipses to a 
future work (C. Beard et al. 2022, in preparation) that employs 
a full photodynamical model (e.g., Carter et al. 2012; Mills & 
Fabrycky 2017). 

After performing these steps, TTVs were detected. We 
phase-folded the individual transits (Figure 4) after taking their 
TTVs into account (Figure 5). Figure 4 shows the phase-folded 
and binned transit light curves of each planet. Without 
accounting for TTVs, the phase-folded transits would have 
appeared V-shaped as opposed to U-shaped, and would have 
led to inaccurate transit parameters. We fit all the planets 
simultaneously with emcee (Foreman-Mackey et al. 2013). 
We initialized 128 walkers near the best-fit model from 
imfit. We ran the MCMC for 50000 links and assessed 
convergence using the Gelman-Rubin potential scale reduction 
factor (Gelman et al. 2014). It dropped to below 1.02, 
indicating good convergence. The resultant posterior distribu- 
tion is summarized in Table 10, while Figure 4 shows the best- 
fit transit models. 

We note that the initial detrending of the light curve and 
isolation of transit windows depends crucially on both a good 
knowledge of the TTVs and the transit durations. We therefore 
iterated the whole process outlined in this section twice to 
ensure convergence. In the Appendix, we present a search for 
additional transiting planets in this system. 


5. Transit-timing Variations 


We modeled the observed TTVs with full N-body integra- 
tions of the orbits. As we will show below, at least some of the 
planets of TOL1136 are likely locked in mean-motion 
resonances. In this case, the TTV signal cannot be adequately 
described by the combination of well-known analytic formulae 
for the near-resonant (Lithwick et al. 2012) and individual- 
conjunction (chopping; Deck & Agol 2015) TTVs based on 
perturbation theory. The TTVs in a fully resonant system 
showing nonlinear dynamics cannot be treated in the same way 
(Agol et al. 2005; Nesvorny & Vokrouhlicky 2016). 
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We integrated the orbits using a symplectic integrator (Deck 
et al. 2014; Wisdom & Holman 1991) with a constant time step 
of 0.1 days, considering only the Newtonian gravitational 
interactions between the six planets and the central star all 
treated as point masses. Over the observational baseline of a 
few years, any relativistic precession should be negligible and 
hence ignored. The model transit times were computed as 
described in Fabrycky (2010) by finding the minima of the sky- 
projected star-planet distances. During this iteration for finding 
transit times, the system was integrated using a fourth-order 
Hermite integrator (Kokubo & Makino 2004). The system 
was initialized using values for the planet-to-star mass ratio, 
orbital period P, eccentricity e, argument of pericenter w, and 
time 7, of inferior conjunction nearest to the epoch 
BJD = 2458680, which was converted to the time of pericenter 
passage т via 2m(T. — т) /Р = Eo — esinEg with Eo = 


1-е (5 
2 arctan tan 
1--е 4 


and the longitudes of ascending nodes were held fixed at т/2 
and 0, respectively. The mass ratios and osculating orbital 
elements were converted to Jacobi coordinates using the 
interior mass in Kepler’s third law as in Rein & Tamayo (2015; 
see their Section 23) and the Wisdom (2006) correction for 
the difference between real and mapping Hamiltonian was 
applied once at the beginning of integration as in Deck et al. 
(2014). The sky plane was chosen to be the reference plane, 
with respect to which the arguments of pericenters and the line 
of nodes were defined. The ascending node was defined with 
respect to the 4-Z-axis chosen to point toward the observer. The 
transit-timing code was implemented in JAX (Bradbury et al. 
2018) to enable automatic differentiation with respect to the 
input parameters (see also Agol et al. 2021) and is available 
through GitHub.” 

The N-body transit time model т(0) as described above was 
used to sample from the posterior probability distribution for 
the model parameters 0 conditioned on the observed transit 
times D = {t;}, p(0|D) сс p(D|0)p(0). We adopted the following 
log-likelihood function: 


Б» Б - mer | егеді, (1) 
2 о 


і 


=) . The orbital inclinations 


Inp(D|0) = 


i 


which is based on the assumption that the observed transit 
times are drawn from the independent and identical Gaussian 
distributions around the model values, with variances c? 
estimated from the modeling of transit and Rossiter-McLaugh- 
lin data in Sections 3 and 4 (Table 8). The residuals of transit 
time fitting did not show clear evidence for any non- 
Gaussianity in the tails of the distributions, as has been seen 
in some other works (Jontof-Hutter et al. 2016; Agol et al. 
2021). 


5! We note that this conversion is different from what is adopted in the 
TTVFast code (Deck et al. 2014), which performs the conversion following 
the Hamiltonian splitting defined by Wisdom & Holman (1991). The difference 
comes from the arbitrariness of how to split the motion into nonperturbed (i.e., 
Keplerian) and perturbed parts, and the resulting mappings between the 
coordinates and orbital elements differ slightly by an amount on the order of 
magnitude of the planet-to-star mass ratio. This difference is well below the 
stated uncertainties of any of the parameters, but it matters when one tries to 
reproduce the TTV signal. 

52 https: / [github.com/kemasuda/jnkepler; in this work, we used commit 
6caclc2. 
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Figure 4. TESS light curve phase-folded and binned after removing the measured TTV in TOI-1136. The red curves are our best-fit transit models. Simultaneous 
transits (where two planets transit the host star) were removed before making this plot. 
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Table 2 

Priors adopted in the TTV Modeling 
Parameter Prior 
Planet/Star Mass Ratio (0,5 x 107^) 
Orbital Period (days) ЖР, — 0.5, Py + 0.5) 
Orbital Eccentricity (0, 0.4) 
Argument of Pericenter UO, 27) 
Time of First Inferior Conjunction (days) ИТ — 0.1, To + 0.1) 


Note. Wa, b) is the uniform distribution between a and b. Symbols Ро and To 
denote the linear ephemeris computed from the observed transit times for each 
planet. The argument of pericenter was wrapped at 27. 


We adopted a prior probability distribution function p(0) 
separable for each model parameter, as summarized in Table 2. 
The sampling was performed using Hamiltonian Monte Carlo 
and the No-U-Turn Sampler (Duane et al. 1987; Betancourt 
2017) as implemented in NumPyro (Bingham et al. 2018; 
Phan et al. 2019). We ran four chains in parallel until we 
obtained at least 50 effective samples for each parameter, and 
the resulting chains had the Gelman-Rubin statistic of 
R « 1.05 (Gelman et al. 2014). 

During the TTV posterior sampling, we did not impose any 
requirement for long-term dynamical stability. Instead, we 
imposed a stability requirement in post-processing, as will be 
described in the next section. The planetary parameters 
reported in Table 10 will be based on the stable TTV posterior 
samples. 


6. Dynamical Modeling 
6.1. Stability Analyses 


After examining the posterior distribution of our TTV 
analysis, we realized that many of the posterior samples would 
experience orbital instability on relatively short timescales. As 
TOI-1136 is about 700 Myr old, it should be stable on similar 
timescales. However, with the TTV data in hand, the TTV 
analysis alone may not be able to pin down the system's 
configuration (with >30 parameters) to the island of stability 
that the real system resides. Near MMR the system is 
dynamically rich, and a small change in the system parameters 
may lead to a very different dynamical behavior. This is 
especially true considering the fine structure of second-order 
resonance and the relatively short TTV baseline of the current 
TESS data. 

We therefore proceeded to trim down the posterior samples 
by removing TTV solutions that quickly become unstable. We 
employed the Python package REBOUND (Rein & Liu 2012). 
We used the built-in nercurius integrator, which is a hybrid 
integrator similar to Mercury by Chambers (1999). mer- 
curius makes use of the symplectic Wisdom-Holman 
integrator WHFast (Wisdom & Holman 1991) when planets 
are far away from each other, and switches to the high-order 
integrator IAS15 (Rein & Spiegel 2015) whenever it detects a 
close encounter within a user-defined distance. We switched 
the integrator when any two planets are less than four mutual 
Hill radii from each other. 

We integrated all the posterior samples from Section 5 for 
1 Myr. We acknowledge that this is much shorter than the 
system's age of ~700 Myr. The choice of 1 Myr was a 
compromise between computation time and gauging the long- 
term stability of the TTV solutions. We did not include tidal 
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effects, which may begin to manifest on timescales longer than 
1 Myr. We removed posterior samples that were flagged as 
unstable by REBOUND. Planets in these systems experienced 
collisions or became unbound. 

To quantify the stability of the remaining posterior samples, 
we further examined the orbital architectures after a 1 Myr 
integration. Using the orbital period of the innermost planet b 
as a proxy, we show in Figure 6 that some posterior samples 
underwent substantial changes in the orbital architecture, even 
though the system remained technically stable. In some cases, 
the orbital period of planet b underwent order-of-unity changes 
from its initial value. Moreover, the orbital period ratio 
between the innermost planets Р./Р, moved significantly off 
resonance (Figure 6). These systems later experienced orbital 
instability when we integrated them to 10 Myr. To maximize 
the long-term stability of our posterior samples, we kept only 
posterior samples in which (1) P, changed by <1% from its 
initial value and (2) Р./Р, changed by <2% from its initial 
value of 3:2 MMR after 1 Myr of N-body integration. These 
criteria are shown as the orange box in Figure 6. 

About 48% of the original posterior samples remained after 
the selections just described. All of our subsequent analyses 
were based on this "stable" posterior sample. Table 10 
summarizes this stable posterior distribution and reports the 
osculating Keplerian elements at the time of reference 
BJD — 2458680. We note that the osculating orbital period 
ratios should not be used to predict future transits or gauge the 
depth of resonance in this system. The osculating orbital 
periods suffer from large uncertainty as they vary rapidly after 
close encounters between planets. Instead, we report the orbital 
period ratios by averaging the osculating orbital period of the 
stable solutions over a time interval of 50,000 days (longer than 
the libration periods of the system; see Section 6.2). The period 
ratios are extremely close to their respective resonance, with 
deviations A = ""/^ р of 6919x107 for be, 


2012-097 x10" for cd, 4413x104 for de, 
4.5 ± 1.6 x 10^ for ef, and 8442.9 x 10 ^ for fg. We 
compare this resonant structure to other known planetary 
systems in Section 7.1. 

Given the limited TTV data and measurement uncertainty, 
we most likely have not located the true island of stability that 
is stable for 700 Myr. Resonant interaction involving several 
planets leads to a finely structured and complex dependence of 
the system's dynamical evolution on the initial parameters. A 
small change in the system configuration may lead to a very 
different dynamical behavior. A similar situation was encoun- 
tered by Gillon et al. (2017) in their early analysis of 
TRAPPIST-1. Most of their TTV solutions became unstable 
on very short timescales (~0.5 Myr). Only years later, when 
TTVs were observed over a longer time span, did Agol et al. 
(2021) find solutions for the orbital architecture of TRAPPIST- 
] that are stable for at least 50 Myr. With this in mind, we 
encourage follow-up transit observations of TOI-1136. 

We also tracked which of the TOI-1136 planets were 
dislodged from resonance first. As shown in Figure 7, planets e 
and f (7:5 second-order MMR) seem to be a weak link in the 
resonant chain: they were the first to be removed from 
resonance in more than 68% of the unstable solutions. This is 
theoretically expected because second-order resonant interac- 
tions are weaker than first-order interactions by another factor 
of orbital eccentricity (e E where k is the order of the MMR; 
Murray & Dermott 1999) and have thinner libration widths in 
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Figure 5. Observed TTVs of the planets in TOI-1136, best-fit TTV model (red curve), and 20 dynamically stable posterior samples (blue curves). All data came from 
TESS observations except for the last transit of planet d, which came from our RM measurement. TTVs from neighboring planets are anticorrelated. The superperiods 
are estimated to be 210,000 days, which is much longer than the current observational baseline. Instead, the TTVs are driven by the libration of the resonant angles 


mi + may 2/3 А 
(Nesvorny & Vokrouhlicky 2016). The libration periods were estimated by P; ғ Ra (222) (Agol et al. 2005; Nesvorny & Vokrouhlicky 2016; Goldberg 


m 


et al. 2022) to be between ~700 days and ~5000 days, with the shortest period corresponding to the bc pair. The observed TTV show variations on similar timescales. 


the semimajor axis (see Figure 9). It has also been suggested 
that many second-order resonances formed by convergent disk 
migration may in fact be overstable (Goldreich & Schlichting 
2014; Xu & Lai 2017) and easily disrupted. 


6.2. Generalized Laplace Resonance 


In this section, we investigate whether TOI-1136 planets are 
indeed in MMR rather than being near resonance by chance. 
The hallmark of true MMR is the libration of the relevant 
resonant angles. For a planetary system near resonance, one can 
decompose the Hamiltonian into the Keplerian, resonant, and 


secular terms (Murray & Dermott 1999). The generalized 
coordinate for the resonant interaction is the resonant angle. For 
two-planet systems, the resonant angle ф takes the form: 


фр = 4№ — Pà + Cp — аул). (2) 


where p апа q are positive co-prime integers, and |p — q| is the 
order of the resonance. The mean longitude А is the sum of the 
mean anomaly M, the longitude О of the ascending node, and 
the argument of pericenter w. The angle w is defined as О + w. 
Following D' Alembert's rule, тт 2 can be an integer combina- 
tion of cc, and wz such that the sum of the coefficients is p — q. 
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Figure 6. Fractional change in the orbital period for planet b (left) and the orbital period ratio between planets b and c (right) for our TTV posterior samples after a 1 
Myr N-body integration (Section 6.1). Posterior samples in which the orbital period of planet b moved more than 1% from its initial value are also those that broke 
away from the resonant configuration (P./P, deviated from 3:2 MMR). We removed these systems (outside the orange box) as they quickly became unstable upon 


longer integration. 


f 


Figure 7. Relative fractions of TOI-1136 planets that became dislodged from 
resonance first in our dynamical integration. Dynamical instability often ensues 
after breaking resonance. Planets e and f (the only second-order MMRs in TOI- 
1136) are usually the first to become dislodged due to the weaker strength of 
second-order MMRs. Together, they departed from resonance first in more than 
68% of the posterior samples that became unstable within 1 Myr. 


The strength of the MMR is proportional to е? 4) For a 
system in true two-body MMR, 4» librates around a libration 
center with limited amplitude, as opposed to circulating 
between 0 to 27. 

Several combinations of wm, and о are allowed by 
D'Alembert's rule, especially for higher-order MMRs (Murray 
& Dermott 1999). Exploring all of them can be cumbersome 
and redundant. Sessin & Ferraz-Mello (1984) suggested a 
canonical transformation such that two-body resonance can be 
described by a single mixed pericenter angle (see also Henrard 
et al. 1986; Wisdom 1986; Batygin & Morbidelli 2013b; 
Hadden 2019): 


f e sina, + g e2 sin w 


(3) 


@\2 = arctan 


, 
f еу сов ол + g e2 COS co» 


where f and g are the coefficients of the disturbing function (see 
the tabulated values in, e.g., Lithwick et al. 2012). Petit et al. 
(2020) used this mixed angle to investigate two-body MMR 
and found it useful for probing the resonant angles in K2-19: a 
system with high eccentricities and limited TTV data. We 
adopt this mixed pericenter angle formulation to analyze the 
two-body resonances in TOI-1136. 

When more than two planets are involved in MMR, one can 
generalize the resonant angle. One can simply subtract the two- 


10 


body resonant angles (Equation (2)) of neighboring pairs to 
remove any dependence on w. For a concrete example, 
consider TOI-1136 b, c, and d: 


Poe = 2Ap — ЗА, + We, (4) 
Фа x А, = 2А4 Tue (5) 
Ped = Pre — Фа = 2Аь — 4А, 2. (6) 


A perceptive reader might point out that the coefficients are 
no longer co-prime and that we should divide by 2. We chose 
not to do so following the suggestion of Siegel & Fabrycky 
(2021). The benefit of keeping the original coefficients is 
that the preferred libration centers for a three-body MMR are 
now near 180° in this formulation. For example, in Kepler-60, 
Gozdziewski et al. (2016) defined the three-body resonant angle 
Poca = Ay — 2А. + Ag. Gozdziewski et al. (2016) reported a 
libration center of ~45°. The underlying two-body MMRs are 
5:4 and 4:3; pca should have been фь.а = 4A, — 8А. + 4۸, in 
the formulation of Siegel & Fabrycky (2021). Correspondingly, 
bca) the libration center, should have been 180°. The 
significance of a libration center of 180° is perhaps best 
understood in the most famous example of Laplace’s resonance 
between the inner three Galilean moons, Io, Europa, and 
Ganymede (e.g. Sinclair 1975). The libration of фес = 
Aj—3Ag--2Ag around 180° ensures that whenever two 
satellites have a close encounter, the third satellite is far away, 
by either 90° or 180°. Such a resonant configuration minimizes 
three-body conjunctions and chaotic interactions, and hence 
enhances the overall stability of the system. This geometric/ 
phase relation holds true even for systems that have experienced 
long-range deviation from two-body orbital period commensur- 
ability (e.g., Kepler-221 Goldberg & Batygin 2021). 

One can extend this process to construct resonant angles 
when more planets are involved. In Table 3, we list the various 
resonant angles for TOI-1136. Before describing the results, we 
highlight an effect that can shift the libration centers. For a 
chain of planets in resonance, their mutual interactions change 
the topology of the Hamiltonian, especially when there is a 
nonadjacent first-order MMR. New libration centers can 
emerge that are shifted away from 180° (e.g., Siegel & 
Fabrycky 2021). A system can be captured in one of the 
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Table 3 
Resonant Angles in Stable TTV Posterior 


Resonant Angle Fraction in Libration 


Libration Center Libration Amplitude’ 


Two-body Resonant Angles 


фь=2Аь — ЗА, + Whe 2 100% 179?1 + 1?5 9°6 + 1°5 
фа=А — 2a + Wea 100% 176°7 + 6°8 14°6 + 6°6 
фа-224 — ЗА, + Bue 100% 180°5 + 1°5 1723 + 7°7 
het=SAve — 75у + 20р 91% 182°1 + 74 36° + 13° 
Фь=2Ау — ЗА + бу 100% 180°3 + 0 19° + 15° 
Three-body Resonant Angles 

фьеа=2Аь — 4X + 2047 99% 196° + 15° 19° + 9° 
Qcae- 1c — 4А + ЗА, 93% 163° + 30° 45° + 22° 
dacr=4Aa — ПА + ТА 51% 173° + 37° 64° + 13° 
Qefe=5Ae — 11A + 6А, 5896 143? + 51° 69° + 19° 
Four-body Resonant Angles 

Фьсае=2Аь — 5А. + 6A4 — ЗА, 76% 24° + 36° 44° + 19° 
Фоа 1А, — 8A, + 14А, — 7А 36% —7° + 39° 72° + 6° 
dacte=4r04 — 16А, + 184 — 6A, 5% 


Note. 1: The libration amplitude is defined as A = XM — (@))? (Millholland et al. 2018; Siegel & Fabrycky 2021). 2: А are the mean longitudes of each planet. 


According to the D' Alembert rule, the longitudes of pericenters w of both planets involved in a mean-motion resonance could contribute to the resonant angles. However, 
with a canonical transformation, the two-body resonance is dependent on just the mixed angle: лә = arctan(fe, sin w; + ge; sin cz) /(fe, cos ол + ge, cos wz) (see 
Section 6.2 for more details). For three-body resonances and above, the lowest-order resonant angles are independent of w. 3: We did not reduce the coefficients to be co- 
prime, following the suggestion by Siegel & Fabrycky (2021); in this way, the three-body resonant angles librate near 180°. 


possible libration centers depending on the order of which 
planets are captured into resonance (Delisle 2017). For 
example, Kepler-223 is in a 3:4:6:8 resonant chain (Mills 
et al. 2016). The bd pair (6:3 = 2:1) and the ce pair (8:4 = 2:1) 
are both examples of nonadjacent first-order MMRs. The three- 
body libration centers were hence shifted to 168? and 130? in 
that system (Siegel & Баһтуску 2021). Delisle (2017) 
suggested that the observed configuration is perhaps most 
consistent with Kepler-223 c and d having been captured into 
an MMR before e and b. Fortunately (or sadly), there is no 
nonadjacent first-order MMR in TOI-1136, so one need not 
worry about (or cannot take advantage of) this effect. 

We integrated the stable TTV solutions from Section 6.1 
forward in time for 50,000 days with REBOUND. We recorded 
the various resonant angles of TOI-1136 listed in Table 3. We 
identified systems in which the resonant angles are clearly 
circulating (¢ varied by much more than 27). Then, to identify 
the librating solutions, we calculated the mean of the resonant 
angles during this 50,000 day period. We also computed the 
libration amplitude using the formula in Siegel & Fabrycky 
(2021) and Millholland et al. (2018): 


2 2 
A= "TOXC - (Фу, 


where ($) is mean of the resonant angle, and М is the 
number of resonant angles sampled. If the libration of the 
resonant angle is sinusoidal in shape and sampled regularly in 
time, then A corresponds to the amplitude of that sinusoid. We 
adopted a generous definition of libration: a system is in 
libration if the amplitude is less than 90?. We can see in Table 3 
that most libration amplitudes are much smaller than this 
threshold. 

Figure 9 summarizes the relationships between the various 
resonant angles and the fraction of librating solutions for each 


(7) 
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angle. We found that the various resonant angles involving 
only first-order resonance have a high probability of libration in 
our stable TTV solutions. The fraction is close to unity for the 
two-body angles, and steadily drops as we move up the 
resonance ladder from two-body resonance to multibody 
resonance. The inner four planets b, c, d, e (Фьсае) һауе а 
76% probability of being a resonant chain. Moreover, the 
libration centers are almost always near 0 or 180? (Table 3) as 
found by Siegel & Fabrycky (2021). The only exceptions are 
the resonant angles involving planets c and d (2:1 MMR). 
Beauge (1994) showed that the topology of the phase space of 
the two-body 2:1, 3:1, n:1 MMR permits two libration centers 
that are shifted from 180? (asymmetric libration; Beauge et al. 
2006). The shifts increase with orbital eccentricity. A planetary 
system may adopt one of these libration centers, or chaotically 
shift between them if the libration amplitude is large enough. 
This was confirmed in our convergent disk migration 
simulations (Section 6.3 and the first panel of Figure 8): 
resonant angles involving TOI-1136c and d are shifted from 
180? by asymmetric libration. 

In contrast to the first-order MMR, the resonant angles that 
involve the only second-order MMR (planets e and f, 7:5) have 
significantly reduced probabilities of libration. Second-order 
MMR, by nature, is much weaker and much more localized in 
the phase space than first-order MMR (see Figure 9 and Murray 
& Dermott 1999). In about 9% of our TTV solutions, the 
second-order resonant angle Qer alternates between circulation 
and libration (Figure 8 lower panel) Alternation between 
libration and circulation is a hallmark of chaos and has been 
previously identified in Kepler-36 (Carter et al. 2012). 
However, we strongly suspect that e and f are indeed in a 
7:5 second-order MMR. In our stable TTV solutions, planets e 
and f do have a ~91% chance of being in two-body libration. 
The observed orbital period ratio differs from 7:5 by only 
4.5 + 1.6 x 10 ^; it seems very unlikely to be coincidental. See 
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Figure 8. Top row: the orbital configurations of TOI-1136 just after our convergent disk migration simulation of TOI-1136 (top left), a TTV solution with apsidal 
antialignment between neighboring planets (top center), and a TTV solution without apsidal antialignment (top right). The dotted lines indicate the pericenters of each 
planet. A classical prediction of convergent disk migration (e.g., Batygin 2015) is that neighboring planets should have antialigned pericenters (except for the 
asymmetric libration of cd in 2:1 MMR; see text). A significant fraction of our TTV solutions conform to this prediction (see Figure 19). The evolution of the two- 
body resonant angles of these solutions librate near 180? over the next 50,000 days (middle panel). However, other TTV solutions are far from apsidal antialignment. 
Planets e and f (7:5 second-order resonance) in these solutions often show chaotic behavior where their two-body resonant angle фе, can oscillate between a state of 


libration and circulation (gray line in the bottom panel). 


Bailey et al. (2022) for a dynamical exploration for the 
observed and expected period ratio of pairs of planets locked in 
second-order MMR. Our current TTV solutions of TOI-1136 
are often chaotic on short timescales, with some Lyapunov 
times on the order of 10? days. Again, we suspect that with the 
current TTV data, we have not located the true island of 
stability in the phase space. The measurement uncertainty is 
particularly obvious for second-order MMR that has a thinner 
libration width in the phase space (right panel of Figure 9). 
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We examined the dominant periodicities of the observed 
TTV. For a near-resonant, circulating system, the TTV occurs 
on the timescale of the “superperiod” Р, = 1/|p/P» — 4/Р,| 
(Lithwick et al. 2012). In contrast, for truly resonant systems, 
the TTV should vary on the timescale of the libration period 


-2/3 
Р, N Pw (= for two-body resonance (Nesvorny & 


Vokrouhlicky 2016; Goldberg et al. 2022). We estimated both 
P, and P, in TOI-1136. As the period ratios are so close to 
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Figure 9. Left: ladder of the resonant angles involving increasingly more planets. We recorded the resonant angles (Table 3) in the stable TTV solutions for 50,000 
days. The fraction of the TTV posterior sample in which the specific resonant angle librates is shown in the bracket. The resonant angles that involve the second-order 
resonance of planets e and f (7:5 MMR) have a significantly reduced probability of libration. Second-order resonances have narrower libration widths compared to 
first-order resonances particularly toward low eccentricities (right panel, calculated with Equation 8.76 in Murray & Dermott (1999)). Our TTV analyses, with 
measurement uncertainty, likely has not located the solution to the true island of stability that real system resides. 


ratios of small integers (Section 6.1), the superperiods P, are 
typically longer than 10° days for TOI-1136. On the other 
hand, the estimated libration periods Р, are between about 800 
and 5000 days (from be to fg) based on Equation (2) of 
Goldberg et al. (2022). The existing TTV data clearly show 
variations in the shorter timescales of Р, (Figure 5). For a more 
empirical test, we applied a Lomb-Scargle periodogram to the 
two-body resonant angles фьс to Pr, in our ТТУ posterior 
solutions. P; indeed span a range of 700 to 5000 days. This is 
further evidence that TOI-1136 planets are in resonance rather 
than near resonance. 


6.3. Convergent Disk Migration 


Simulating the formation of resonant-chain planetary 
systems with disk migration can constrain the disk density 
and turbulence, as well as the order of planets that are captured 
into resonances (e.g., Hühn et al. 2021). Previous works (Xu & 
Lai 2017) have shown that it is much more challenging to form 
a second-order MMR than first-order MMR through disk 
migration. If the disk migration were turbulent or simply rapid, 
a planet pair could have easily skipped a second-order 
resonance and become locked in nearby first-order resonances. 
We can leverage this difficulty of forming the observed second- 
order 7:5 MMR of TOI-1136 ef to constrain the properties of 
TOI-1136's protoplanetary disk. 

We experimented on three prescriptions of disk migration for 
TOL1136 (see schematics in Figure 10). Our first set of 
simulations follow the prescription of Cresswell & Nelson 
(2006), Baruteau et al. (2014), Pichierri et al. (2018), and Hühn 
et al. (2021). Type-I migration was applied to all the planets 
simultaneously. The rate of migration on each planet was 
calculated based on the planetary properties and their current 
locations in the protoplanetary disk (for details see Section 3 of 
Pichierri et al. 2018). This procedure was implemented in the 
type I migration routine of REBOUNDx (Tamayo et al. 
2020). Crucially, to halt the migration and prevent planets from 
plunging into the host star, we included an inner edge of the 
disk in the simulations. The existence of an inner edge in a 
protoplanetary disk at the corotation radius is theoretically 
expected (e.g., Ghosh & Lamb 1979; Ostriker & Shu 1995). 
Observationally, the inner edge may also be responsible for the 
decline of sub-Neptune occurrence inward of 10 days (e.g., 
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Terquem & Papaloizou 2007; Lee & Chiang 2015). The 
location of the inner edge was set to be 0.05 au (near the 
current orbit of TOI-1136 b), with a transition region of 0.01 au 
over which the migration torque is smoothly reversed to mimic 
the effect of the pressure bump. Planet b was initialized 5% 
outside its currently observed orbit. The other planets were 
initialized with orbital separations such that each pair has a 
period ratio 2% wider than their currently observed resonances. 
This is to represent in-situ formation of the planets followed by 
short-scale (—0.1 au) migration. The planetary masses were 
taken from the stable TTV posterior samples. The only 
exceptions are planets e and f, which were assigned mass 
ratios of 0.9 < q < 1.1 and the same mass scale from ТТУ 
solutions. As suggested by Xu & Lai (2017), having a mass 
ratio near unity maximizes the chance of establishing and 
maintaining a second-order MMR. The planets had initially 
circular orbits and randomized arguments of pericenter and 
mean anomalies. The main tunable parameters in this 
simulation are the surface density of the protoplanetary disk 
at 1 au (У, а) and the scale height h=H/R. We assumed 
У= Dy au a ‘> and varied У au between 10 and 10* gom ? 
uniformly in logarithmic space. Thus, the simulated disks have 
surface densities that vary between about 1/200 and 10 times 
that of the minimum-mass solar nebula (Weidenschilling 1977; 
Hayashi 1981 > 4,221700 g cm 2)? 3 h was randomly chosen 
between 0.01 and 0.1 and was assumed to be constant 
throughout the disk (no disk flaring). For easier comparison 
with the typical disk lifetime of a Sun-like star (~3 Myr; see, 
e.g., Andrews 2020), we converted [> au, Л] to the decay 
timescale of the semimajor axis and orbital eccentricity (та, Te] 
using the equations in Pichierri et al. (2018). The whole system 
was evolved for 3 7,; visual inspection of the time evolution 
confirmed that all planets have had ample time to complete 
migration and settle into MMR (Figure 11). 

Our second prescription of disk migration is widely used in 
the literature: e.g., Tamayo et al. (2017) employed this method 
to successfully simulate the formation of TRAPPIST-1 (Gillon 
et al. 2017). In this prescription, Type-I migration was only 
applied to the outermost planet. The benefit is that all 
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See also the minimum-mass extrasolar nebulae, whose surface densities 
have substantial variation between different systems (Chiang & Laughlin 2013; 
Dai et al. 2020). 
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Figure 10. Summary of our disk migration simulations (Section 6.3). We experimented with three prescriptions of disk migration: (1) we applied Type-I migration to 
all planets simultaneously with a disk edge (left column). (2) Type-I migration was only applied to the outermost planet (a scenario that may happen in transition disk; 
middle column). (3) Similar to the first case except that the planets migrated from beyond 1 au as opposed to the 0.1 au in the previous two cases (right column). The 
top row shows the schematics for each mode of migration. The second row shows the results of the simulation in terms of migration timescales in 7, and К = т,/т, 
compared with the typical disk lifetime (~3 Myr for Sun-like stars; Andrews 2020). The blue filled symbols are simulations that formed analogs of TOI-1136 where 
planets are in their observed resonances particularly with planets e and f in a second-order 7:5 MMR. The red hollow symbols are systems that have failed (usually e 
and f skipped 7:5 and became locked in a nearby first-order MMR). The third row shows the final orbital period ratio between planets e and f. The fourth row shows 
the depth of MMR produced in A at the end of the simulations. The gray area indicates the observed A between planets b and c. In general, the first prescription, short- 
scale (from 0.1 au) disk migration with a disk edge, is the most robust at producing systems of TOI-1136. It can deposit systems deeper in MMR with A < 10 * as 
was observed in TOI-1136. The migration process could be completed quickly within the typical disk lifetime. 


encounters between the planets are now convergent: the inner 
planets do not migrate until they are captured in resonances 
with the outer planets. This prescription may seem contrived; 
however it may be the case in transition disks (Espaillat et al. 
2014) where the inner gas disk is starting to disperse (see 
schematic Figure 10). There can be a time at which only the 
outermost planet is still embedded in a gas disk and 
experiences Type-I migration. We dynamically evolved the 


system using RI 


EBOUND with the WHFAST integrator (Rein & 
Liu 2012). The effect of Type-I migration was implemented 
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using the modify_orbits_forces routine in REBOUNDx 
(Tamayo et al. 2020). As we are migrating just one planet, we 
directly varied 7, uniformly in logarithmic space between 10* 
and 10’ yr. Instead of varying т, directly, we varied К = T4/Te 
between 10 and 1000 (right panels of Figure 10). K bears 
theoretical significance that will be explained shortly. 

Our third prescription is almost identical to the first 
prescription. We applied Type-I migration to all planets 
simultaneously, and we included an inner disk edge. The only 
difference is that the planets were initially placed further out in 
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Figure 11. Simulated Type-I migration where migration was applied to all planets, and there is an inner edge of the disk at 0.05 au. The panels, respectively, show the 
time evolution of the deviation from MMR A, orbital eccentricities, two-body resonant angles, and three-body resonant angles. The inner disk edge halts the migration 
of the planets and turns initially divergent encounters into convergent encounters (first panel). As shown here, the system was quickly captured into a resonant chain 
including the second-order resonance for planets e and f on a timescale of 10^ yr. Once in resonance, eccentricities are excited by resonant interaction, while the 


resonant angles start to librate. 


the disk (>1 au). This prescription specifically investigate the 
ex-situ formation of the TOI-1136 planets followed by large- 
scale migration. 

Figure 11 shows the time evolution of the period ratios, 
orbital eccentricities, and resonant angles in a successful disk 
simulation using the first prescription. Figure 12 displays the 
resultant constraint on the disk properties. The planets were 
locked into their observed MMR on 10s kyr timescales. Once 
in resonance, the resonant angle changed from a state of 
circulation to libration. Even though the planets started on 
circular orbits, resonant interaction can pump up the eccen- 
tricity. Another well-known result of convergent disk migration 


is that the deviation from MMR A = E — ] and the 
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equilibrium orbital eccentricity e are inversely related (e.g., 
Ramos et al. 2017). The inverse relation is determined by the 
ratio between the semimajor axis and eccentricity damping 
timescales К =т„/т„. During disk migration, e is damped 
down by the disk and is pumped up by resonant interaction. 
The equilibrium eccentricity is given by the balance of the e 
pumping and e damping (Terquem & Papaloizou 2019). In a 
Sessin-type resonant Hamiltonian (безіп 4 Ferraz- 
Mello 1984), if we ignore secular interaction and work in the 
limit of small e, the argument of pericenter precesses at a rate 
ш x 1/e for a planet in MMR (see also Laune et al. 2022). For 
a pair of planets to remain in MMR, the period ratio has to 
deviate away from MMR (A increases) such that the 
conjunctions shift spatially in pace with the precession of 
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Figure 12. Properties of the protoplanetary disk that formed TOI-1136. Left: total disk surface density and scale height (Л =H/R) of our disk migration simulation. 


The successful simulations (blue solid symbols) suggest that TOI-1136 likely had a lower total surface density (Utotai1 au 


< 1000g cm ?) than the minimum-mass 


~N 


solar nebula (MMSN; Hayashi 1981). The slower migration in a lower density disk facilitated the capture into resonance particularly the 7:5 second-order resonance. 
Right: MMSN of TOI-1136 constructed directly from the TTV-measured masses using the method in Dai et al. (2020). TOI-1136 falls close to the Kepler multiplanet 
systems with an estimated solid surface density of Хат au 50 g ст 2. The two panels together suggest an enhanced dust-to-gas ratio of at least >0.05 within the 
innermost | au possibly due to radial drift and gas disk dispersal (e.g., Birnstiel et al. 2010; Gorti et al. 2015; Cridland et al. 2016). 


pericenters: фуз = qn, — рп, + c ~ 0. Our disk migration 
simulations recovered this general behavior (bottom row of 
Figure 10). A smaller K= 7;/7,, slower damping of orbital 
eccentricity, leads to larger equilibrium e and smaller deviation 
from MMR A (Figures 10 and 14). To reproduce the observed 
A of ^10 (gray area in Figure 10), К--т./т, has to be 
smaller than about 100. 

We carried out about 200 simulations for each prescription. 
This was not an exact number as some realizations became 
unstable. The results are summarized in Figure 10. We consider 
a simulation successful if all six planets get locked into their 
observed MMR with no more than 0.1% deviation; the 
respective two-body and three-body resonant angles are all 
librating. The most common failure mode is that planets e and f 
skip the weaker second-order 7:5 MMR and become locked in 
the nearby stronger first-order MMR (4:3 and 3:2; see third row 
of Figure 10). Our simulations disfavored the third prescription: 
a long-scale (from 1 to 0.05 au) Type-I migration. None of the 
200 simulations with this prescription managed to form a 
system like TOI-1136. Xu & Lai (2017) found that the capture 
into a second-order resonance is more likely with slower 
migration (see their Equation (44)). There is a paradox here that 
if the planets experience long-scale migration, their migration 
rate must be high enough so that they can arrive at the observed 
0.05 au separation before the disk dissipates after ~3 Myr. On 
the other hand, the weak 7:5 second-order resonance is easily 
skipped during fast migration. Even though in some realiza- 
tions planets e and f get initially captured into 7:5 MMR, 1 au 
to 0.05 au is such a long journey that the perturbations that 
form the other planets eventually disrupt the weak 7:5 MMR. 

Our two short-scale (0.1 au) migration prescriptions both 
abundantly produce TOI-1136 analogs (Figure 10). However 
the second prescription, migrating only the outermost planet, 
seems less likely. To form analogs of TOI-1136, the second 
prescription often requires slower migration with timescales of 
several Myr that often exceed the typical disk lifetime (second 
row of Figure 10). One may argue that in transition disks, the 
gas surface density is low enough that Type-I migration is also 
significantly slower. However, a transition disk is short-lived 
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leaving little time for the migration to deposit the planets deep 
in resonance. In particular, the innermost planets have to wait 
for the outer planets to be captured into resonance sequentially 
before resonant interaction starts acting on it. Our simulations 
very rarely deposit planets b and c to the observed 1074 level 
from perfect resonance (bottom row of Figure 10). 

Our first prescription, a short-scale (0.1 au) Type-I migration 
on all planets with a disk edge, seems to be the more likely 
scenario. As shown in Figure 10, the first prescription can 
produce systems like TOI-1136 (including 7:5 MMR) even 
with a rapid Type-I migration of 7, = 10*106 yr. This is 
thanks to the inner edge of the protoplanetary disk, which 
slows down and even reverses the effective migration (Kretke 
& Lin 2012; Masset et al. 2006). The disk edge stalls the inner 
planets at the edge and thereby allows planets further out to 
catch up and join the resonant chain (Izidoro et al. 2017). As 
shown in the top panel of Figure 11, even though some planet 
pairs initially underwent divergent migration, all planet pairs 
eventually switched to convergent migration and got locked 
into MMR. Moreover, as all planets migrated simultaneously 
and captured into resonance quickly, they are deposited deeper 
in resonance after the simulation (A can be as low as 1075; 
bottom row in Figure 10). Such deep resonances better match 
the observed TOI-1136 system. Within the limitations of the 
Type-I migration prescriptions of Cresswell & Nelson (2006), 
Baruteau et al. (2014), and Pichierri et al. (2018), our 
successful disk migration simulation translates to a proto- 
planetary disk no denser than ~1000 g cm ? at lau 
(Figure 12). This is comparable but lower than the surface 
density of the MMSN (241700 с cm >; Hayashi 1981). 

Another signpost of convergent disk migration is that 
neighboring planets in MMR have antialigned arguments of 
pericenters (e.g., Batygin & Morbidelli 2013a). This is a robust 
prediction of convergent disk migration as it does not depend 
on the initial conditions. Antialigned pericenters have been 
observed in some of the known resonant chains (e.g., 
TRAPPIST-1; Agol et al. 2021). For TOI-1136, our disk 
migration simulations produced antialigned pericenters for 
neighboring planets (see the top left panel of Figure 8). A 
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significant fraction of our TTV posterior samples is indeed 
consistent with an antialigned configuration (top center panel of 
Figure 8). However other TTV solutions are not apsidally 
antialigned (top right panel of Figure 8). On a population level 
(Figure 19), our TTV solutions are suggestive of apsidal 
antialignment; however more TTV data are needed to confirm 
this trend. 

Macdonald & Dawson (2018) proposed that post-formation 
eccentricity damping alone could also produce a resonant chain 
of planets. In our limited exploration of this possibility, we 
could only deposit the inner two or three planets of TOI-1136 
into resonance. The other planets, which have much longer 
tidal timescales (see Section 6.4), showed negligible evolution 
within a 700 Myr age. We argue that post-formation 
eccentricity damping may explain pairs or triplets of resonant 
planets; however it struggles to explain a six-planet resonant 
chain such as TOI-1136. Some other process, e.g., Type-I 
migration, is required to initialize the planets close to 
resonance. The observed orbital architecture of TOI-1136, 
particularly the depth of MMR and the 7:5 second-order MMR, 
is most consistent with the scenario of a short-scale (0.1 au), 
Type-I migration with an inner disk edge. 


6.4. Resonant Repulsion 


After the protoplanetary disk dissipates, a resonant chain of 
planets in a system like TOI-1136 may experience planetesimal 
scattering (e.g., Chatterjee & Ford 2015), orbital instabilities 
followed by giant impact collisions (e.g., Izidoro et al. 2017; 
Goldberg & Batygin 2022), secular chaos (e.g., Petrovich et al. 
2019), and tidal dissipation (e.g., Lithwick & Wu 2012), all of 
which could move the system off resonance. If the system is 
lucky, it may evade giant impacts and planetesimal scattering; 
however some amount of tidal resonant repulsion (Papaloizou 
б Terquem 2010; Lithwick & Wu 2012; Batygin & 
Morbidelli 2013a; Delisle & Laskar 2014; Pichierri et al. 
2019) seems unavoidable. Resonant repulsion is well under- 
stood for a pair of planets in first-order MMR: tidal damping of 
both planets’ eccentricities causes A to rise, taking the system 
further from MMR. This is not to be confused with a simple 
divergence of orbits due to the faster tidal orbital decay of the 
inner planet. In resonant repulsion, the outer planet moves 
outward. The underlying physics is almost identical to the e— 
A-K relationship we described in Section 6.3. Again, when the 
resonant interaction dominates and in the limit of small e, the 
argument of pericenter precesses at a rate w ос 1/e. To stay in 
MMR, the period ratios of two resonant planets must positively 
deviate away from MMR to catch up with the ever faster 
precession of the pericenter. In short, as e gets damped by tides, 
w precesses faster, and A has to increase to maintain the MMR. 
This process can continue until the resonance is broken. Again 
Kepler-221 is a great example (Goldberg & Batygin 2021). 

Most of the Kepler multiplanet systems are not near first- 
order MMR. There is only a small overabundance just wide of 
first-order resonances and a lack of planets just short of them 
(Fabrycky et al. 2014). See also Figure 15. A number of works 
have explored whether this preponderance of wide-of-reso- 
nance systems could be produced by resonant repulsion (e.g., 
Lee et al. 2013; Silburt & Rein 2015). The general conclusion 
is that, with only eccentricity tides, resonant repulsion is too 
slow to explain the entire Kepler sample. Millholland & 
Laughlin (2019) pointed out that obliquity tides may solve this 
problem by enhancing the rate of tidal dissipation. Regardless 
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of the source of dissipation, the long-term asymptotic behavior 
of resonant repulsion is the same, as long as the process does 
not break the resonance (or the Cassini state for the case of 
obliquity tides; Batygin & Morbidelli 2013a). The long-term 
asymptotic behavior obeys a power-law relation: A ос (t/7,)!/3 
(e.g., Equation (26) of Lithwick & Wu 2012). 

We simulated the resonant repulsion for TOI-1136. Тһе 
initial conditions are our disk migration simulations that 
successfully locked all six planets of TOL1136 into their 
observed MMR (Section 6.3). We integrated these systems 
forward in time in REBOUNDx (Tamayo et al. 2020). We used 
the symplectic WHFAST integrator (Wisdom & Holman 1991). 
We included tidal damping on all planets using the mod- 
ify orbits forces routine in REBOUNDx. We parame- 
terized 7, using the equilibrium-tide expression (Murray & 


Dermott 1999) 
2 m Р 
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where и is the mean motion of the planet; О is the tidal quality 
factor; kz is the tidal Love number; т, is the stellar mass; and 
ть, ть, and a are the planetary mass, radius, and semimajor 
axis, respectively. 

To guide our discussion, we first examine the theoretical 
behavior of resonant repulsion for each pair of planets in TOI- 
1136 (Lithwick & Wu 2012). This calculation assumes that the 
planets are only in pairwise first-order resonance. According to 
their Equation (26), A grows as (t/7;)!/? with a proportionality 
that changes with the planetary parameters and the relevant 
resonance. The process of resonant repulsion is independent of 
the absolute scale of т, as long as the system is maintained in 
resonance. This is why Goldberg & Batygin (2021) were able 
to use a т, of just 10 yr to speed up their numerical 
investigation of Kepler-221. The situation is more complicated 
for a resonant chain of planets: the effect of tidal damping on 
individual planets will be transmitted to other planets by their 
resonant interaction. Resonant repulsion could proceed for all 
resonantly locked planets even though tides may only operate 
strongly on the inner, larger planets. 

With six planets in TOI-1136, each of which may have a 
different reduced tidal quality factor Q' = Q/k», there are too 
many possibilities to consider. For simplicity, we assumed that 
all planets have the same Q' but different 7, given by 
Equation (8). In this case, planets b, c, and d have comparable 
T, ~5 Gyr if the planets have Neptune-like О’ ~ 3 x 10^ (e.g., 
Banfield & Murray 1992; Zhang & Hamilton 2008). Planets e, 
f, and g have т, values that are longer by at least 2 orders of 
magnitude. However, we also tried to decrease the Q' of planet 
b by 2 orders of magnitude than the other planets. This 
possibility was entertained because planet b is plausibly rocky 
(1.9R5), which would make it much more dissipative than a 
gaseous planet. We also explored the possibility that planets d 
and f may have Q' smaller than the other planets by 2 orders of 
magnitude. The motivation is that d and f are the largest 
planets; perhaps their radii are inflated by the heat of obliquity 
tidal dissipation. 

We integrated TOI-1136 for about 100 7, of the most 
dissipative planet to determine the asymptotic behavior of 
resonant repulsion. The qualitative behavior is similar regard- 
less of the choice of Q' and is shown in Figure 13. The 
theoretical A ос (1/т.)!/3 behavior held up well even though 
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Figure 13. Time evolution of the deviation from MMR A in our resonant repulsion simulations of TOI-1136 (Section 6.4). We dynamically evolved the resonant 
chains generated from our convergent disk migration simulations (Section 6.3) after including tidal dissipation. We note a characteristic behavior in which A grows 
with (1/7;)!/? is seen for all planet pairs as long as they remain in a resonant chain. The x-axis is plotted with т, of the most dissipative planet. Even though tidal 
dissipation may be concentrated on the most dissipative planet (usually planet b), resonant repulsion occurs on all planet pairs as long as they remain a resonant chain. 


Depending on the masses of the planets, the rate of resonant repulsion is typically on the order of 107° (solid lines). However, the observed deviations are on 


(/ re3 
the order of 10~* (horizontal dotted lines). The intersection between solid lines and horizontal lines tells us the number of т, cycles that have elapsed іп the 700 Myr 
lifetime of TOI-1136. As shown in the plot, the intersection happened at small (t/7,)!/3 indicating minimum tidal evolution since formation. We described in the text 
that we rule out a few scenarios that would enhance the rate of tidal dissipation. 


TOI-1136 is in a resonant chain rather than a resonant pair, for Table 4 

which the theory was originally derived (Batygin & Rate of Resonant Repulsion 

Morbidelli 2013a; Lithwick & Wu 2012). We experimented Simulated A 

with т, values between 10” and 10" yr, and confirmed that the т, of after 700 Myr Largest 
qualitative behavior stayed the same. We concluded using Scenario Planet b ! Observed A 
Equation (26) of Lithwick & Wu (2012) that planet pairs bc, Earth-like (Q/io x 100) 12My 14% <008% 
cd, de, and fg would deviate from MMR A by about 0.0005 to Mars-like (О/К» ~ 1000) 120 Myr 0.7% <0.08% 
0.004 after 1 7e. A would double these amounts after 87, given Neptune-like (О/ ~ 30,000) 3.8 Gyr 0.1% 50.08% 
the 1/3 power law. REBOUNDx simulations revealed consistent 

rates of resonant repulsion of 0.0006 to 0.004 for various planet Note. 1. Deviation from MMR A after 700 Myr of resonant repulsion. 
pairs (Figure 13). The precise values depend on the TTV- Reported here is the planet pair that shows the fastest deviation. See the text for 
measured masses. details. 


The observed deviations from MMR (A) coupled with an 
age estimate for the system can be used to put constraints on 
the tidal quality factor Q’ of the planets (Brasser et al. 2022; 
Lithwick & Wu 2012). Compared to the other known systems 
with resonant chains, TOI-1136 is young, with an estimated 
age of 700 Myr. In Figure 13, we plotted the evolution of A as 
a function of time; again note the the long-term asymptotic 
behavior is A œ (t/7,)!/3. In theory, the intersection of the 
currently observed A (horizontal dashed lines) and the resonant 
repulsion A « (t/ 7,) 5 power law (solid lines) could provide 
an empirical estimate of 7.. However, the А œ (t/7,)!? 


convergent migration, the system was deep in resonance (A 
can be as low as 10 7 in our disk migration simulations) with 
large orbital eccentricities (e ~ 0.1). As tides operate, eccen- 
tricies get damped and resonant repulsion drives the system 
toward larger A. We plotted the measured e and A constraints 
from our TTV analyses in Figure 14. The relatively high e and 
small A in our TTV solutions are consistent with the very end 
of disk migration or the very start of resonant repulsion. In 
other words, TOI-1136 has undergone very minimal resonant 
repulsion and still records the orbital architecture from disk 


relation only holds asymptotically for long-term evolution migration. Even the most dissipative planet in TOI-1136 likely 
(Figure 13). Due to other terms in the Hamiltonian, the early- has 7, that is at least 700 Myr if not much longer. For example, 
time behavior deviates significantly from a perfect (т/т„)!/5 if all of the planets in ТОІ-1136 have Neptune-like 
relation. Nonetheless, we can see that intersection happened Q' x 3 x 104, т, would be at least 4 Gyr, and one would not 
early on with (1/7,)? < 1. In other words, the 700 Myr old expect to see significant resonant repulsion in its 700 Myr 
TOI-1136 has barely undergone a single 7, of the most lifetime. In contrast, most Kepler near-resonant TTV systems, 
dissipative planet. We can hence rule out an Earth-like or Mars- typically a few Gyr old, have A ~ 1% and damped eccentricity 


like О’ for planet b (1.9Rẹ). If b had a terrestrial Q’ of 1000 ex 0.02 (e.g., Hadden & Lithwick 2014). They have likely 
(Murray & Dermott 1999), 7, would have been ~120 Myr. undergone many cycles of tidal damping thanks to perhaps 
About 5 т, cycles have elapsed in TOI-1136’s lifetime, A obliquity tides (Millholland & Laughlin 2019) or other 


would have been significantly higher than the observed value. mechanisms. 

We summarize a few representative cases in Table 4. Planetesimal scattering can also induce deviations from 
The constraints on the orbital eccentricity from our TTV MMR (Chatterjee & Ford 2015). One can put an upper limit on 

analysis also shed light on the progress of resonant repulsion. the integrated amount of planetesimal scattering based on the 

In the e-A space (Figure 14), each planet follows an evolution extremely deep resonances observed in TOI-1136. For Kepler- 

track that is anticorrelated in e and A. The underlying physics 223, Moore et al. (2013) found that there could not have been 

was explained at the start of this section. Soon after the more than one Mars mass of orbit-crossing planetesimals or the 
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Figure 14. The evolution of deviation from MMR A and the orbital eccentricity provide a different perspective on resonant repulsion (see Figure 13). A resonant 
planetary system starts with substantial eccentricity and is deep in resonance right after disk migration (upper left corner of this A—e parameter space). Over time the 
orbital eccentricities are damped by tides, and the deviation A from MMR increases. Qualitatively, this is because the precession of pericenter scales inversely with 
eccentricity in a Sessin-type Hamiltonian (Sessin & Ferraz-Mello 1984). To maintain resonance, the orbital period has to deviate from the perfect integer ratio to catch 
up with the precession. The measured constraints on A and e from our TTV analysis are shown by the error bars. The high-e and low-A TOI-1136 system likely has 
not undergone much resonant repulsion since formation. Brasser et al. (2022) made a similar plot for TRAPPIST-1. With e < 0.005 and A > 1%, TRAPPIST-1 is a 
mature (a few Gyr old) system that has likely evolved to the bottom right of this A—e parameter space. 


systems would have been pulled out of resonance. Similarly, 
Raymond et al. (2022) investigated the same question for 
TRAPPIST-1. TOL1136 may be amenable to a similar 
investigation, which is left for a future work. 


7. Discussion 
7.1. A System Deep in Resonance 


TOI-1136 is a deeply resonant planetary system. We now 
compare it with other known multiplanet systems. We only 
included planets discovered by the transit method in this 
comparison because orbital periods are much more precisely 
measured in transit surveys than in other types of surveys. We 
did not include the TESS objects of interest (TOI; i.e., Guerrero 
et al. 2021) because TOIs typically have much shorter 
observational baselines; hence the orbital periods are not as 
precisely measured as in the Kepler mission. Moreover, many 
TOIs have not been confirmed yet. 

Figure 15 shows that TOI-1136 stands out as one of dozen 
planetary systems with orbital periods extremely close to 
MMR. Near-resonant Kepler multiplanets typically deviate 
positively from MMR by about 1 to 2% (Fabrycky et al. 2014). 
However, the planets orbiting TOI-1136 have A that are 
roughly 2 orders of magnitude smaller according to our 
analyses (Section 6.1). The other planetary systems with 
similarly low A are also resonant-chain systems, such as 
Kepler-60 (GoZdziewski et al. 2016; Jontof-Hutter et al. 2016) 
and Kepler-223 (Mills et al. 2016). 

Another metric for identifying resonance was proposed by 
Goldberg & Batygin (2021): B = pn, — (p + q)no + qn, = 
Dig It quantifies how fast the three-body resonant angle 
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changes with time. This metric is useful for picking out 
planetary systems that are in generalized Laplace resonance, in 
which case the resonant angle librates and B is small in 
magnitude. In TOI-1136, the values of B for the neighboring 
triplets bcd, cde, def, and efg are all smaller than in the general 
Kepler multiplanet systems by at least 1 order of magnitude 
(Figure 16). Again, the planetary systems with similarly small 
B are those with resonant chains: Kepler-221 (Goldberg & 
Batygin 2021), Kepler-223 (Mills et al. 2016), Kepler-60 
(Gozdziewski et al. 2016), etc. Even without a TTV analysis, 
the depths of resonance among the six TOI-1136 planets seem 
extremely unlikely unless there is some underlying physical 
process that drove the planets into resonance. 

Our TTV and dynamical analyses (Section 4 and 6.2) 
provided further evidence for a resonant chain in TOI-1136. 
We showed that the planets of TOI-1136 display TTV on 
timescales that are more consistent with the libration period of 
the resonant angles (700—5000 days) than the superperiod or 
the circulation of the resonant angles (710,000 days). More- 
over, our stable TTV solutions predominantly showed the 
libration of the various resonant angles (Figure 8, 9 and 
Table 3) with libration centers near the theoretically predicted 
values (Siegel & Fabrycky 2021). 


7.2. Planet ef (7:5 MMR) is the Weakest Link 


The ef pair is the only second-order resonance in TOI-1136. 
Previous investigation has shown that second-order MMR is 
both much more difficult to form from disk migration and more 
easily disrupted than first-order MMR (Xu & Lai 2017). This is 
because a second-order MMR, compared to first-order 
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Figure 16. B is another metric for identifying multiplanet systems in generalized Laplace resonance (Goldberg & Batygin 2021). В = [рт — (р + g)n2 + qna 


, where 


п = 2т/Р is the mean motion, and p and 4 are co-prime integers. For resonant systems, the resonant angle is librating; hence its time derivative B should be small in 
magnitude. Plotted here is B normalized by the average mean motion <n> of the Kepler multiplanet systems (blue), TOI-1136 (red), and other known resonant chains 
(orange). B values in TOI-1136 are similarly low as the other resonant-chain systems. TOI-1136 is by far the most observable resonant-chain system with a V-band 


magnitude of 9.5. 


resonance, is suppressed by a factor of the orbital eccentricity e. 
Moreover, the width of the second-order MMR in the phase 
space is much thinner (Figure 9 of Murray 4 Dermott 1999). 
Mah (2018) showed that planets b, c, and d of TRAPPIST-1 
(second- and third-order resonance) were often the first to be 
displaced from resonance in their dynamical simulations. 
Similarly, our dynamical modeling of TOI-1136 (Section 6.1) 
indicated that planets e and f are often the first to be dislodged 
from resonance. Dynamical instability often ensues after the ef 
pair is removed from the resonant chain. 

We further experimented with the possibility that the ef pair 
is in the nearby 3:2 and 4:3 first-order MMR despite a period 
ratio that is close to 7:5 commensurability. One notable 
example is Kepler-221 (Goldberg & Batygin 2021); the planets 
are in a Laplace resonance even though their pairwise period 
ratios deviated by >10% from small integer ratio. For TOI- 
1136, we analyzed the resonant angles Фе, baer, and Perg in 100 
random draws of the stable TTV solutions. In all of these 
solutions, the resonant angles фер, Oder, and Pefg are circulating 
when computed with 3:2 or 4:3 MMR. A 7:5 second-order 
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resonance for planets e and f is the simpler and preferred 
solution. 

We also explored the possibility that there is an additional 
planet between planets e and f such that the planets are in a 
chain of 5:6:7 first-order MMR. The existence of such a planet 
would eliminate the need that planets e and f are in the much 
weaker 7:5 second-order MMR. In Appendix B, we show that 
both a systematic transit search and a careful visual inspection 
were not able to detect this hypothetical planet. Moreover, we 
calculated the mutual Hill radius for the 5:6:7 configuration. 
The planets are separated by only six mutual Hill radii even if 
the hypothetical middle planet is only about ІМ. Such tight 
packing is seen in <0.5% of all Kepler multiplanet systems and 
may compromise the overall stability of the system (Pu & 
Wu 2015). Furthermore, including this hypothetical planet did 
not lead to an improved TTV solution. All of these results are 
against the possibility of another planet between planets e 
and f. 

Therefore planets e and f are likely indeed in a 7:5 second- 
order MMR. This represents a weak link in the resonant chain 
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of TOI-1136 and may threaten the overall stability as the 700 
Myr old system continues to mature. In the Kepler multiplanet 
sample, there is an overabundance of planets just outside first- 
order resonance (Fabrycky et al. 2014). However there is no 
noticeable feature near second-order resonance (except perhaps 
5:3; Steffen & Hwang 2015). The discovery of TOI-1136 
shows that second-order MMR can be produced in at least 
some protoplanetary disks, as suggested by Xu & Lai (2017). If 
so, does it mean that the observed paucity of second-order 
MMR in mature planetary systems is due to the dynamical 
fragility of such a configuration? Izidoro et al. (2017) was 
puzzled that, in order to reproduce the observed fraction of 
resonant systems in Kepler, at least 75% (or even 95%) of their 
simulated, initially resonant planetary systems must become 
unstable. However, only 5096-6096 of their first-order MMR 
chains became unstable. Maybe the inclusion of the weaker 
second-order MMR could increase the rate of orbital instability. 
We note that in the revised models of Izidoro et al. (2021), the 
fraction of unstable planetary systems could reach 95%. 
Moreover some of their simulated planetary systems contained 
second-order MMRs. 


7.3. Comparison with Other Resonant Chains 


TOI-1136 joins a handful of known planetary systems with 
a resonant chain: GJ 876 (Rivera et al. 2010; Mullholland 
et al. 2018), TRAPPIST-1 (Gillon et al. 2017; Luger et al. 
2017; Wang et al. 2017; Agol et al. 2021), TOI-178 (Leleu 
et al. 2021), Kepler-80 (MacDonald et al. 2016), Kepler-60 
(Gozdziewski et al. 2016), K2-138 (Christiansen et al. 2018), 
Kepler-223 (Mills et al. 2016), and Kepler-221 (Goldberg & 
Batygin 2021). K2-72 (Crossfield et al. 2016), V1298 Tau 
(David et al. 2019) as well as the Kepler systems labeled in 
Figure 15 might also have resonant chains, pending further 
analysis. Tejada Arevalo et al. (2022) argued that V1298 Tau 
cannot be in resonance based on stability considerations. 

TOI-1136 is the second known resonant-chain system with a 
well-established age as young as a few hundred million years. 
The other system is Kepler-221, with an age of about 600 Myr 
(Goldberg & Batygin 2021). The rest of the resonant-chain 
systems are at least several Gyr old or have no precise age 
estimates. TOI-1136 and Kepler-221 seem to have had 
disparate evolution tracks despite similar ages. In Kepler-221, 
although the pairwise orbital period ratios (1.765 and 1.829) 
are farther from commensurability than in TOI-1136, the three- 
body resonant angle changes so slowly (small B; Figure 16) 
that the resonant angle is most likely librating. The interpreta- 
tion offered by Goldberg & Batygin (2021) is that Kepler-221 
underwent rapid tidal resonant repulsion, possibly with the help 
of obliquity tides. Goldberg & Batygin (2021) estimated that a 
total of 70007, must have elapsed for the system to reach the 
current state of 1096 off resonance. On the other hand, TOI- 
1136 has barely moved from perfect orbital period commen- 
surability. One possible explanation is that the conditions for 
capturing planets into a Cassini state (Millholland & Laughlin 
2019) were simply not available for TOI-1136. Its resonant 
repulsion has to proceed with the much slower eccentricity 
tides. We will return to this point in Section 7.6. Based on the 
preceding argument, Kepler-223 (not to be confused with 
Kepler-221; Mills et al. 2016) may represent the future of TOI- 
1136. Kepler-223 is about 6 Gyr old, and its four transiting 
planets are likely in a four-body resonant chain that only 
involves first-order MMRs. Despite its 6 Gyr age, Kepler-223 
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seems to have avoided giant impact collisions, resonant 
repulsion, and planetesimal scattering, any of which could 
have induced deviations from MMR. Its orbital period ratios 
are still deep in resonance (1.3336, 1.5015, and 1.3339). 

TOI-1136 also has the first known resonant chain that has a 
second-order MMR between neighboring first-order MMRs. 
Kepler-29 b and c have a period ratio that deviates from a 9:7 
MMR at a 10 * level (Fabrycky et al. 2012; Jontof-Hutter et al. 
2016); however existing TTVs could not determine if the 
system is in resonance nor its dynamical origin (Migaszewski 
et al. 2017; Vissapragada et al. 2020). TOI-178 b is near a 
second-order 5:3 MMR with planet c (Leleu et al. 2021). 
However, the period ratio is shorter than expected if the system 
was resonant (1.95 day versus 1.91 day). Leleu et al. (2021) 
suspected that tidal dissipation might have broken planet b 
away from resonance. In TRAPPIST-1, it is also the case that 
the inner three planets are close to third-order (8:5) and second- 
order (5:3) MMR; Agol et al. (2021) showed that the three- 
body resonant angle involving b, c, and d is likely librating. 
However, they could not tell if the two-body resonant angles 
were also librating. Nonetheless, the presence of a 490 day 
superperiod in the TTV of TRAPPIST-1 suggests the 
circulation of the two-body resonant angle. One may be 
tempted to suggest that these innermost planets were initially in 
first-order MMR, but later disrupted by a crossing of the disk 
edge or tidal evolution (Huang & Ormel 2022). TOI-1136 is a 
very rare case—possibly unique among the known systems— 
of a resonant chain with a second-order MMR between first- 
order resonances. 


7.4. The Disk that Formed TOI-1136 


The second-order resonance between TOI-1136 e and f 
allows us to place stringent constraints on TOI-1136's 
formation environment. Planets e and f most likely started 
with an initial period ratio close to 1.4 such that they did not get 
captured into the nearby, much stronger 3:2 first-order MMR. 
Xu & Lai (2017) showed that the successful capture and 
stability of a second-order MMR is facilitated by lower initial 
orbital eccentricity, a planet mass ratio т»/т\ close to unity, 
and most importantly slower disk migration. 

In our disk migration simulations (Section 6.3), the disks that 
successfully locked e and f into a 7:5 MMR all had a lower 
total surface density (hence lower migration rate) compared to 
ће MMSN (514510006 ст 2). Іп comparison, Hiihn et al. 
(2021) used a very similar prescription of disk migration with 
an inner disk edge to constrain the formation of Kepler-223, 
which only contains first-order MMR (Mills et al. 2016). Hiihn 
et al. (2021) noted that Kepler-223 could form from convergent 
disk migration with a wider range of disk properties: the disk 
surface density can be a few times denser than the MMSN but 
still lock all planets of Kepler-223 into a resonant chain. 

The rate of disk migration allowed us to constrain the total 
disk surface density of TOI-1136's protoplanetary disk. Our 
analyses in Section 6.3 suggested that the TOI-1136 planets 
formed mostly in situ followed by short-scale migration. If so, 
one can also constrain the solid surface density by spreading 
out the masses of the planets into their local feeding zones. We 
computed the minimum-mass extrasolar nebula of TOI-1136 
using the TTV masses following the method in Dai et al. 
(2020). TOI-1136 joined the other Kepler multiplanet systems 
with a similar solid surface density of Хұоңа, 12422206 cm ? 
(Figure 12). The total surface density X, au < 1000 g cm ? and 
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the solid surface density together suggest an enhanced dust-to- 
gas ratio of 20.05 within the innermost 1 au of TOI-1136. This 
is higher than the typically assumed value of 0.01 in the 
interstellar medium, and may suggest radial drift of dust and 
early gas disk dispersal (e.g., Birnstiel et al. 2010; Gorti et al. 
2015; Cridland et al. 2016). 

Previous disk migration simulations placed meaningful 
constraints on disk turbulence (Adams et al. 2008; Rein & 
Papaloizou 2009; Hühn et al. 2021). Turbulence may increase 
the libration amplitudes of the planets captured in resonance 
and even disrupt the resonance if the turbulence is strong 
enough. We did not include turbulence in our convergent disk 
migration simulations in Section 6.3, because the libration 
amplitudes of TOI-1136 are still poorly constrained by the 
available TTV data. Recent work by Jensen & Millholland 
(2022) indicated that typical methods for inferring libration 
amplitudes of resonant planetary systems can be strongly 
biased by measurement uncertainties. We defer a discussion of 
the libration amplitudes to a future work where the libration 
amplitudes are better constrained. 

TOI-1136 has a highly coplanar planetary system that is also 
well-aligned with the host star. The fact that all six (potentially 
seven) planets transit already hints at a low mutual inclination. 
According to our transit modeling (Section 4), assuming all 
planets transit the same hemisphere of the host star, the 
measured orbital inclination implies a mutual inclination of 
1?1. The planet with the most discrepant orbital inclination is 
planet b at 86.447021^. The other five planets all have orbital 
inclinations around 8925, If we assume that the planets have the 
same longitudes of ascending node, which can be tested with 
future transit duration variation analyses, the dispersion of 
orbital inclinations (0°15) is proxy for their mutual inclination. 
Previous works have also found that the innermost planet of a 
multiplanet system often has the largest orbital inclination, 
likely due to an equipartitioning of the angular momentum 
deficit (Steffen & Coughlin 2016; Dai et al. 2018; Petrovich 
et al. 2019; Weiss et al. 2018b). For comparison, the 
TRAPPIST-1 planets have even lower mutual inclinations of 
about 0°04 (Agol et al. 2021). 

Our RM measurement of TOI-1136 revealed a planetary 
system that is well-aligned with the rotation of the host star. 
The sky-projected stellar obliquity А is 5° + 5°, and the stellar 
obliquity is less than 28° with a 95% credible level. Hirano 
et al. (2020) measured the stellar obliquity of TRAPPIST-1, 
and they also found evidence for a well-aligned planetary 
system. Spalding & Batygin (2016) pointed out that planets 
may couple with the oblateness of their host star, the 
differential nodal precession may induce a mutual inclination 
of ~2W. Hence, for TOI-1136, the measured low mutual 
inclination and low stellar obliquity corroborate each other: if 
the stellar obliquity were high, a large mutual inclination would 
have been generated by the differential precession. Batygin 
(2015) suggested that if a protoplanetary disk has substantial 
axial asymmetry, capturing planets into MMR during disk 
migration is much more difficult. If TOI-1136 had a stellar 
companion, a stellar flyby event (e.g., Xiang-Gruess 2016) or 
the perturbation from the companion (Batygin 2012) may 
induce disk warp, axial asymmetry, and primordial misalign- 
ment that are detrimental to the formation of a resonant chain 
like TOI-1136. No spectroscopic, visual, blended, or comoving 
stellar companions were found for TOI-1136 (Section 2). We 
suggest that TOI-1136, which is 700 Myr old, still preserves 
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the pristine orbital architecture formed by a slow migration in 
an isolated disk with no primordial misalignment or disk 
asymmetry. 


7.5. Mass and Radius 


Figure 17 shows the masses and radii of the TOI-1136 
planets along with the theoretical mass-radius relationships 
from Zeng et al. (2016) and Chen & Rogers (2016). We also 
plot archival mass measurements from the NASA Exoplanet 
Archive.” Using the model by Chen & Rogers (2016), which 
takes the age, mass, composition, and insolation level of a 
planet into account, the required mass of the H/He envelope 
increases from 0.1% for the smallest planet b up to about 15% 
in mass for the largest planet d. 

TOI-1136 is about 700 Myr old. The innermost planet b 
should have experienced extensive photoevaporation for 
hundreds of Myr (Fulton et al. 2017; Owen & Wu 2017; Zhang 
et al. 2022). With a core mass of about ЗМ, planet b does not 
have a deep enough gravitational potential well to prevent 
photoevaporation (see the self-consistent hydrodynamic simu- 
lations of photoevaporation by Wang & Dai (2018)). On the 
other hand, the more massive and more distant planets in TOI- 
1136 should experience sequentially weaker photoevaporation. 
A dynamically quiet, multiplanet system like TOI-1136 is a 
good test bed for photoevaporation theory. Given the delicate 
orbital architecture, the orbital distances likely did not change 
significantly since formation. There have not been any giant 
impact collisions that could have removed the gaseous 
envelopes (Inamdar & Schlichting 2016). The planets are 
subject to the same XUV spectrum other than scaling with their 
orbital distances. By comparing the extent of the mass loss for 
each planet, TOI-1136 offers an opportunity to probe the 
variance of the efficiency of photoevaporation (Owen & 
Campos Estrada 2020). Future observations of outflowing 
material via Lya or metastable He observations (e.g., Spake 
et al. 2018; Zhang et al. 2022) coupled with hydrodynamic 
modeling (e.g., Wang & Dai 2021) may shed light on this 
issue. 

In Figure 17, we color-coded mass measurements from TTV 
and RV analyses separately. TOI-1136 conforms to the 
previously suggested trend that TTV planets tend to have 
lower masses than RV planets of the same radii (e.g., 
Steffen 2016). One possible explanation is that for near- 
resonant systems (majority of the TTV sample), strong e—mass 
degeneracy may bias the TTV measurements toward lower 
values (Lithwick et al. 2012). However, for TOI-1136, a 
resonant TTV case, such a degeneracy is minimal (Nesvorny & 
Vokrouhlicky 2016). Instead, both the amplitude and periodi- 
city of resonant TTV contain information on the masses of the 
planets independently. Hence the е-таѕѕ degeneracy does not 
seem to be a convincing explanation for TTV planets’ lower 
densities. Alternatively, TTV planets might have originated 
from beyond the water snow line where conditions are more 
conducive for accreting a thick atmosphere (e.g. Lee & 
Chiang 2016). Another possibility is obliquity tides: Millhol- 
land & Laughlin (2019) argued that obliquity tides might 
inflate the radii of near-resonant systems as the tidal dissipation 
goes into heating the planets. We discuss the obliquity tides in 
more detail in the next section. 


ч https: //exoplanetarchive.ipac.caltech.edu 
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Figure 17. Measured masses and radii of the TOI-1136 planets and other 
planets on the NASA Exoplanet Archive. We distinguish the mass 
measurements from RV (blue) and TTV analyses (gray). It was pointed out 
that TTV planets may have systematically lower masses than RV planets of the 
same radii (e.g., Steffen 2016; Hadden & Lithwick 2017; Mills & Mazeh 2017). 
Although there is still a large mass uncertainty based on the existing TTV data 
set, TOI-1136 planets tend to have lower densities than the RV planets even 
though the e-mass degeneracy does not affect this in-resonant system 
(Nesvorny & Vokrouhlicky 2016). The theoretical mass-radius relationships 
are from Zeng et al. (2016) and Chen & Rogers (2016). TOI-1136 b, being the 
innermost, lowest-mass planet, seems to have lost substantial H/He after 
700 Myr of photoevaporation. The outer, more massive planets are generally 
consistent with having ~2% to 15% their mass in H/He. 
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7.6. Obliquity tides inflating d and f? 


Strong tidal dissipation due to obliquity tides may offer 
enough tidal damping to explain both the observed resonant 
repulsion in Kepler multiplanet systems (Millholland & 
Laughlin 2019) and a possible radius inflation of near-resonant 
planets (Millholland 2019). Obliquity tides can be much more 
dissipative than eccentricity tides. When a planet maintains a 
nonzero obliquity (Cassini State 2 being the most favorable; 
Colombo 1966; Peale 1969), the tidal bulge will move in the 
corotating frame and lead to significant dissipation. The capture 
of a system into a Cassini state (secular spin—orbit resonance) 
both excites and sustains a nonzero planetary obliquity. This 
resonance happens when the nodal precession frequency 
matches the spin precession frequency. Resonant planetary 
systems are more likely to be in a Cassini state. This is because 
during their convergent disk migration, the nodal precession 
frequency sweeps through a range of frequencies as the 
semimajor axes of the planets change, allowing a crossing of 
the frequencies. 

TOI-1136 d and f have larger radii (>4Rẹ) than the other 
planets (<3R,), to some extent discrepant with the previously 
noted trend of intrasystem uniformity of multiplanet systems 
(Millholland et al. 2017; Wang 2017; Weiss et al. 2018a). 
Could their larger radii be caused by obliquity tides? As 
Millholland & Laughlin (2019) and Millholland (2019) 
suggested, obliquity tides can lead to both resonant repulsion 
and radius inflation of resonant planets. The TOI-1136 planets 
are currently still deep in resonance with deviations from MMR 
A on the order of 10 ^. The system has undergone minimal 
resonant repulsion since formation (Section 6.4). Even 
assuming the planets have moved by a A of 10 2, the 
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corresponding tidal heating luminosity spread in the system's 
age of 700 Myr is about 10/6 W. This only accounts for 10 ^ of 
the bolometric insolation planet d receives from the star 
(102° үу), According to the Modules for Experiments in Stellar 
Astrophysics simulation by Millholland (2019), this level of 
additional heating due to tides could not inflate the planetary 
radius by more than 10%. This is not sufficient to inflate the 
radii of d and f to >4Rç if they were initially similar to the 
other planets with radii «3R4. The radii of d and f may require 
another explanation such as slower photoevaporation, or dusty 
outflows (Wang & Dai 2019; Gao & Zhang 2020), which is 
testable with a near-infrared transmission spectrum from the 
James Webb Space Telescope (Gardner et al. 2006). In this 
hypothesis, the radius inflation of planets d and f is temporary. 
As the system matures, the radii may drop down to conform to 
the intrasystem uniformity of mature multiplanet systems. 


7.7. A Precursor of Kepler Multiplanet Systems? 


Finally, we place TOI-1136 in the broader context of the 
formation and dynamical evolution of close-in, sub-Neptune 
planets. Figure 18 shows our understanding of where the field 
stands and how TOI-1136 fits in. Planet embryos grow in 
protoplanetary disks. The rate of Type-I migration is propor- 
tional to the masses of the cores (Kley & Nelson 2012). 
Therefore, depending on the rate of core growth, Type-I 
migration may or may not play a significant role in the 
formation of close-in sub-Neptune planets. For systems where 
planetary cores assembled quickly, Type-I migration may 
routinely generate a chain of resonant planets parked at the 
inner edge of the disk. TOI-1136 is an example of this scenario: 
itis a 700 Myr old adolescent planetary system that still records 
the deeply resonant configuration from disk migration. On the 
other hand, in systems with slower core growth, planet 
embryos undergo limited migration and are generally non- 
resonant when the disk dissipates. Otherwise, if the disk is 
turbulent or axially asymmetric (Adams et al. 2008; Batygin 
2012), one would also expect a nonresonant configuration. 
Post-disk assembly of nonresonant planetary systems will 
likely remain nonresonant (right column of Figure 18). 

Fast-forwarding to a ~5 Gyr old mature planetary system, 
TOI-1136 may remain deeply resonant if it only experiences 
negligible dynamical evolution such as orbital instability 
(Izidoro et al. 2017) and planetesimal scattering (Chatterjee 
& Ford 2015). The deeply resonant, 6 Gyr old Kepler-223 
(Mills et al. 2016) may be the future of TOI-1136. So far there 
are about ~10 resonant-chain systems. If TOI-1136 continues 
to undergo mild resonant repulsion and planetesimal scattering, 
it may join the population of near-resonant (A = 196—246), 
multiplanet Kepler systems that show circulating TTV. There 
are about 100-200 such systems that have been discovered by 
Kepler (e.g., Jontof-Hutter et al. 2016; Hadden & Lithwick 
2017). If the future dynamical evolution of TOI-1136 is more 
violent, orbital instability and giant impact collision (Izidoro 
et al. 2017; Goldberg & Batygin 2022) may totally disrupt the 
resonance in TOI-1136. It will end up as a nonresonant 
planetary systems that dominates the mature Kepler multiplanet 
sample. 


8. Summary 


Disk migration may be a common stage of planet formation 
(Ward 1997; Kley & Nelson 2012). If so, many close-in, tightly 
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Figure 18. Schematic showing how TOI-1136 fits into the broader picture of planet formation. It provides an adolescent planetary system that still records the initial 
condition from convergent disk migration. Depending on whether its future dynamical evolution is quiescent or violent, it may stay as a resonant chain like Kepler-223 
(bottom left; Mills et al. 2016), mildly evolve into a near-resonant system (bottom center), or has its resonant structure violently disrupted and become a nonresonant 
Kepler multiplanet system (bottom right). 
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packed, multiplanet systems as observed by Kepler should host 
planets in a chain of mean-motion resonances (e.g., Lee & 
Peale 2002; Cresswell & Nelson 2008; Kley & Nelson 2012). 
In reality, only a small subset of Kepler multiplanets are 
observed near resonance with a typical deviation A of 1% to 
2% (Fabrycky et al. 2014). Over billions of years of dynamical 
evolution, a combination of effects including resonant repul- 
sion (Lithwick & Wu 2012), dynamical instability (Goldberg 
et al. 2022), secular chaos (Wu & Lithwick 2011), and 
planetesimal scattering (Chatterjee & Ford 2015) could lead to 
the slow deviation or the disruption of the migration-induced 
resonance. 

To complete this picture, we present TOI-1136, a young 
planetary system with a resonant chain of six planets. The 
system is so deep in resonance that it probably still preserves a 
"pristine" orbital architecture from convergent disk migration. 
It may be a precursor of many of the Kepler near-resonant 
multiplanets before dynamical evolution eventually dislodged 
the planets from perfect resonance over the Gyr timescale. Our 
observations and dynamical modeling revealed the following 
characteristics of TOI-1136: 


1. TOI-1136 is about 700 Myr old based on gyrochronol- 
ogy, activity indicators, and Li absorption. 

2. A Rossiter-McLaughlin measurement of planet d 
revealed a planetary system well-aligned with the host’s 
rotation with a sky-projected stellar obliquity of 5° + 5°. 
All six planets transit, which implies a low mutual 
inclination between the planets: 1°1 or just 0°15 after 
excluding the most inclined planet, b. 

3. No spectroscopic, adaptive optics, visual, or comoving 
stellar companion was detected for TOI-1136. The low 
stellar obliquity, coupled with the coplanarity, and 
dynamical fragility of a resonant chain of planets, point 
to the formation of TOI-1136 in an isolated disk with no 
stellar flyby, disk warp, or significant axial asymmetry. 

4. There are six transiting planets with each neighboring 
pair showing anticorrelated TTVs. The TTVs are most 
likely driven by the libration of resonant angles (libration 
periods) rather than by the circulation of resonant angles 
(superperiods). 

5. Our TTV analysis revealed the masses of the planets. The 
mass and radius of the innermost and lightest planet, b, 
suggests only a 0.1%-by-mass H/He envelope. This is 
consistent with the expectation of 700Myr of 
photoevaporation. 

6. The orbital period ratios are extremely close to the ratios 
of small integers, with a deviation A — E — 1 оп 


the order of 10 ^. 

7. The closeness to MMR and the libration of the various 
resonant angles suggest that TOI-1136 planets are in 
resonance rather than near resonance. 

8. Planets e and f are close to a 7:5 second-order MMR. 
TOI-1136 is the first known resonant chain with a 
second-order MMR between first-order ММК. Тһе 
weaker and more delicate second-order MMR is much 
more difficult to form in disk migration and more easily 
dislodged from resonance later on. 

9. Our disk migration simulations favor Type-I migration 
with an inner disk edge for TOI-1136. The edge helps to 
halt the migration of the planets and converts divergent 
encounters into convergent ones. To lock the ef pair into 
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a 7:5 second-order MMR, the disk has to be less dense 
than than the MMSN with >, a, < 1000 g ст. 

Our resonant repulsion simulations indicate that TOI- 
1136 has undergone minimum tidal dissipation since its 
formation. Strong tidal dissipation due to a rocky planet b 
or obliquity tides on planets d and f seems unlikely. 


10. 


We encourage additional photometric follow-up observation of 
this system using space-based and ground-based facilities in the 
next few years to refine the dynamical constraints on this 
system. TOI-1136 is also amenable to metastable helium 
observation and transmission spectroscopy that will help better 
understand this young planetary system. 
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Appendix A 
TFOP Observations 


TOI-1136 received a number of follow-up observations 
from the TESS Follow-up Observing Program (TFOP). We 
refer the readers to the full list of observations on ExoFOP. 
We briefly summarize them here. As part of the standard 
process for validating transiting exoplanets and assessing the 
possible systematic errors in the planetary radius due to light 
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from bound or unbound companions (Ciardi et al. 2015), 
TOI-1136 was observed with higher-resolution instruments 
including near-infrared adaptive optics (AO) imaging at 
Palomar, Gemini-North, and Lick, optical speckle imaging at 
Gemini-North Scott et al. (2021), and lucky imaging on the 
AstraLux instrument (Hormuth et al. 2008) at the Calar Alto 
Observatory. The optical observations generally provided 
higher resolution than the NIR observations, while the NIR 
AO generally provided better sensitivity (especially to low- 
mass stars). The combination of the observations in multiple 
filters enables better characterization for any companions that 
may be detected. 

Two reconnaissance spectra were obtained on UT 2019 
December 3 and UT 2020 January 28 with the Tillinghast 
Reflector Echelle Spectrograph (TRES; Fűrész 2008) located on 
the 1.5 m telescope at the Fred Lawrence Whipple Observatory 
(FLWO) in Arizona. TRES is an echelle spectrograph that 
operates in the wavelength range 390—910 nm and has a resolving 
power of 44,000. The TRES spectra were extracted using 
procedures outlined in Buchhave et al. (2010). The TRES spectra 
were also visually inspected. No signs of a composite spectrum 
(blended binary) were found. The TRES spectra were also used to 
derive stellar parameters using the Stellar Parameter Classification 
(SPC; Buchhave et al. 2012) tool. The resultant stellar parameters 
agreed well with our HIRES results presented in Section 2. SPC 
gave Ter = 5775 + 50 К, log g = 4.47 + 0.10, [m/H] = —0.02 + 
0.08, vsini, = 6.7 + 0.5 kms". 

The KeplerCam on the 1.2 m telescope at the FLWO was 
used to catch a transit of planet c on UT 2020 January 25 in the 
Sloan-z band. AstroImageJ (Collins et al. 2017) was used to 
perform aperture photometry and model the predicted event. 
Unfortunately, this observation was performed before the team 
realized there is TTV in TOI-1136. We did not detect the transit 
event. 


Appendix B 
Search for Additional Planets 


We systematically search for a seventh transiting planet in the 
TESS light curve. A BLS analysis did not detect another 
significant signal beyond the six known planets. We performed a 
visual inspection of the TESS light curve, which revealed a 
possible seventh planet in this system. We saw a single transit-like 
event centered at BJD-2457000 — 2435.10 (Figure 21) with a 
duration of about 7.4 hr. Assuming the planets are all on circular 
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orbits, such a transit duration would imply an orbital period of 2.1 
times that of planet g (Р, = Б,(Т, IEY x 2.1). However, after a 
thorough visual inspection, we could not identify another transit 
event of similar depth and duration in the existing TESS light 
curve. We analyzed this planet simultaneously with other planets 
in TOI-1136 following the procedure in Section 4. The transit 
depth implied a planetary radius near 2.5R,, although the data 
could accommodate a planet up to 5А. if the planet is on a 
grazing orbit. We tried to add this seventh planet to our TTV 
model (Section 5). We assumed that TOI-1136.07 followed the 
trend of the other planets and had an orbital period exactly twice 
that of planet g. However, adding this planet to our TTV model 
does not lead to an improvement of the fit. The fit looked almost 
identical visually, and there is no improvement in the Bayesian 
information criterion (Schwarz 1978). We did not include this 
planet candidate in our final analysis. 


Appendix C 
HARPS-N Rossiter-McLaughlin Measurement 


We observed a spectroscopic transit of TOI-1136d using 
HARPS-N (Cosentino et al. 2012; Mayor et al. 2003) mounted on 
the 3.58 m Telescopio Nazionale Galileo (TNG) located on Roque 
de los Muchachos, La Palma, Spain. We observed a transit of 
TOL1136 d on the night starting on UTC 2021 May 14 with 
observations between 21:30 UT and 04:00 UT. The exposure time 
was set to 900s and with an overhead of roughly 20s, the 
sampling was approximately 920s. Due to varying weather 
conditions the signal-to-noise ratio (in order 49) ranged from 
around 80 (in the beginning of the night) to around 40—50. 

We used our HARPS-N transit data to get an independent 
measure for the projected obliquity of TOI-1136 d. We 
sampled the posteriors using MCMC sampling using the 
emcee (Foreman-Mackey et al. 2013) package with the code 
by Hirano et al. (2011) to model the RM effect. We imposed 
Gaussian priors to Ry, a /R,, and i according to the values in 
Table 10. Gaussian priors were imposed to the macro- and 
microturbulence with values of 3.13 + 1 kms" (Doyle et al. 
2014) and 1.04+ 1 kms ! (Bruntt et al. 2010), respectively. 
We let the mid-transit time and the systemic velocity, the 
v sini,, and the sky-projected obliquity A to vary freely. The 
posterior distribution indicates А of 6729 consistent with the 
HIRES measurement. The mid-transit time of this event was at 
BJD of 2459349.525 + 0.005. 
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Figure 19. Eccentricity vectors of adjacent planet pairs in TOI-1136. The colored contours represent a successively higher posterior probability for each planet. The 
gray contours show the relative eccentricity vector between neighboring planets (е; соѕил у — e;cosw;). The relative eccentricity vectors better capture any 
covariance. A classical prediction of convergent disk migration scenario is that adjacent planets should be apsidally antialigned (see Section 6.3 for details). For 
apsidally antialigned solutions, the gray contours tend to be pushed away from the origin. Existing constraints on TOI-1136 may hint at apsidally antialigned 
configurations; however more data are required to confirm this trend. 
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Figure 20. Rossiter-McLaughlin measurement of TOI-1136 d during a transit near ВЛ) = 2459349.525 with HARPS-N. We measured a stellar obliquity А of 6*38°, 


which is consistent with the higher-S/N HIRES measurement. 
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Figure 21. Top: A possible single transit of a seventh planet in TOI-1136 was identified by visual inspection near BJD —2457000 = 2435.10 (left). We are unable to 
confirm this planet; no similarly shaped transit event was seen in the rest of the TESS light curve. Bottom: best-fit transit model of the detrended and binned light 
curve. The nominal transit depth suggests a planetary radius of about 2.5Rẹ. However, many posterior samples are also consistent with a larger planet (584) on а 
grazing orbit. The transit duration is about 7.4 hr. The planets were on a circular, edge-on orbit. The implied orbital period is roughly twice the period of planet g 
(Р, = Р,(Т, /T,)° & 2.1) following the resonant pattern. 
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Figure 22. TESS light curves of TOI-1136 across different sectors (sectors 14, 15, 21, and 22). The quasiperiodic flux modulation due to stellar rotation is clearly 
visible. We removed these variations by fitting a cubic spline (orange curves) to the out-of-transit fluxes (green). The top panels show the detrended light curve and 
mid-transit times of planets if there were no TTVs. 
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Figure 23. Same as Figure 22 for sectors 41 and 48. 
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Table 5 
Keck/HIRES Radial Velocities during a TOI-1136 d Transit near BJD = 2459650.0310 
Time (BJD) RV (ms) RV Unc. (ms) 
2459649.865577 0.02 1.26 
2459649.870924 —2.86 1.32 
2459649.876399 0.34 1.44 
2459649882382 —1.91 1.51 
2459649.888204 —2.63 1.49 
2459649.894153 0.31 1.28 
2459649.90016 —4.02 1.35 
2459649.905831 2.21 1.16 
2459649.91 1004 —1.68 1.24 
2459649.91597 1.32 1.19 
2459649.92105 —1.3 1.23 
2459649.926027 —0.95 12 
2459649.930888 —1.85 1.34 
2459649.935772 0.87 1.24 
2459649.940657 2.97 1.26 
2459649.945529 0.4 1.37 
2459649.950413 0.14 1.31 
2459649.955355 8.23 1.27 
2459649.960807 9.05 1.26 
2459649.966779 8.41 1.18 
2459649.971964 9.3 1.3 
2459649.976848 8.87 1.21 
2459649.98 164 11.07 1.27 
Table 6 
Keck/HIRES Radial Velocities during a TOI-1136 d Transit near BJD = 2459650.0310 Continued 
Time (BJD) RV (ш!) RV Unc. (ms) 
2459649 .986466 8.68 1.19 
2459649.99 1408 10.88 1.18 
2459649.996327 9.78 1.24 
2459650.001685 10.73 1.29 
2459650.007449 7.21 1.22 
2459650.013653 —0.09 1.21 
2459650.019509 —2.09 1.36 
2459650.025967 2.91 1.3 
2459650.032345 0.92 1.31 
2459650.038467 3:15 1.19 
2459650.044567 -5.95 1.35 
2459650.050666 —5.54 1.2 
2459650.056233 —8.68 1.25 
2459650.061568 —10.1 1.24 
2459650.067182 —10.54 12 
2459650.073432 —6.54 1.29 
2459650.079751 —5.96 1.26 
2459650.085874 —12.66 1.25 
2459650.092135 —10.78 1.28 
2459650.098547 —10.01 1.33 
2459650.104531 —8.63 1.2 
2459650.110318 —3.7 1.25 
2459650.116023 —0.24 1.26 
2459650.121822 —2.95 1.21 
2459650.127343 —0.08 1.32 
2459650.132852 —0.92 1.31 
2459650.138465 —2.98 1.28 
2459650.144055 —3.64 12 
2459650.150016 1.72 1.29 
2459650.156208 —0.39 1.22 


3l 


THE ASTRONOMICAL JOURNAL, 165:33 (37pp), 2023 February Dai et al. 


Table 7 
HARPS-N Radial Velocities during a TOI-1136 d Transit near BJD = 2459349.525 
Time (BJD) RV (ms) RV Unc. (ms) 
2459349.400949960109 7421.3 2.2 
2459349.411447130144 7423.6 2.0 
2459349.422129469924 7423.7 2.0 
2459349.433135880157 7422.1 2.3 
2459349.443552029785 7426.2 2.5 
2459349.454326970037 7430.5 32 
2459349.465472249780 7431.9 S. 
2459349.475332879927 7432.9 4.9 
2459349.486998970155 7432.9 4.2 
2459349.497172090225 7429.4 4.3 
2459349.507634540088 7430.7 4.9 
2459349.517553030048 7423.5 5.4 
2459349.529068669770 7420.5 5.4 
2459349.540619030129 7411.2 4.4 
2459349.551579140127 7416.7 3.6 
2459349.561578650028 7420.3 3.3 
2459349.5728 16519998 7419.7 3.2 
2459349.582538269926 7416.3 32 
2459349.593209039886 7420.8 3.3 
2459349.604898279998 7415.0 72 
2459349.615210279822 7426.33 7.5 
2459349.626124090049 7434.2 6.1 
2459349.636806440074 7415.4 5.4 
2459349.647245740052 7430.1 4.5 
2459349.6574767 19934 7436.7 5.6 
2459349.669582610019 7429.7 6.8 
Table 8 
Measured Mid-transit Times of TOI-1136 Planets 

Planet Epoch Mid-transit Times (BJD-2457000) Unc. (days) 
b 0 1684.2689 0.0128 
b 1688.4659 0.0153 
b 2 1692.6029 0.0085 
b 4 1700.9705 0.0193 
b 5 1705.1523 0.0162 
b 6 1709.3189 0.0118 
b 7 1713.4489 0.0213 
b 8 1717.6700 0.0124 
b 9 1721.8301 0.0203 
b 10 1726.0011 0.0109 
b 11 1730.1877 0.0170 
b 12 1734.3482 0.0108 
b 45 1872.0104 0.0072 
b 46 1876.2113 0.0148 
b 47 1880.3842 0.0182 
b 49 1888.7064 0.0152 
b 50 1892.8671 0.0162 
b 51 1897.0779 0.0112 
b 52 1901.2685 0.0096 
b 53 1905.4083 0.0182 
b 56 1917.9092 0.0198 
b 57 1922.1214 0.0164 
b 58 1926.2296 0.0173 
b 177 2422.6729 0.0102 
b 178 2426.8482 0.0080 
b 179 2431.0311 0.0086 
b 180 2435.1835 0.0089 
b 181 2439.3783 0.0135 
b 182 2443.5572 0.0112 
b 222 2610.4300 0.0108 
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Table 8 
(Continued) 
Planet Epoch Mid-transit Times (BJD-2457000) Unc. (days) 
b 224 2618.8061 0.0100 
b 226 2627.1275 0.0166 
b 227 2631.3006 0.0106 
b 228 2635.4951 0.0137 
Note. 1. From HARPS RM measurement (not included in our TTV modeling). 2. From HIRES RM measurement. 
Table 9 
Measured Mid-transit Times of TOI-1136 Planets Continued 

Planet Epoch Mid-transit Times (BJD-2457000) Unc. (days) 
c 0 1688.7211 0.0036 
c 1 1694.9699 0.0023 
c 2 1701.2284 0.0028 
c 3 1707.4861 0.0022 
c 4 1713.7520 0.0032 
c 5 1720.0034 0.0018 
c 6 1726.2605 0.0080 
c 7 1732.5187 0.0022 
c 30 1876.4569 0.0019 
c 33 1895.2336 0.0022 
c 34 1901.4896 0.0029 
c 35 1907.7538 0.0024 
с 37 1920.2739 0.0087 
с 117 2420.9969 0.0018 
с 118 2427.2600 0.0026 
с 120 2439.7757 0.0017 
с 121 2446.0284 0.0018 
c 148 2614.9957 0.0024 
c 149 2621.2509 0.0033 
c 150 2627.5170 0.0022 
c 151 2633.7730 0.0039 
d 0 1686.0671 0.0012 
d 1 1698.5858 0.0011 
d 3 1723.6219 0.0013 
d 4 1736.1413 0.0013 
d 15 1873.8428 0.0012 
d 16 1886.3601 0.0011 
d 18 1911.3936 0.0012 
d 19 1923.9092 0.0012 
d 53 2349.525 0.005! 
d 59 2424.6430 0.0010 
d 60 2437.1649 0.0012 
d 74 2612.4673 0.0012 
d 77 2650.0310 0.0018? 
e 0 1697.7758 0.0022 
e 1 1716.5624 0.0099 
e 2 1735.3536 0.0097 
e 10 1885.7918 0.0044 
e 11 1904.5934 0.0070 
e 12 1923.4102 0.0112 
e 39 2430.9549 0.0034 
e 49 2618.8132 0.0044 
f 0 1699.3854 0.0018 
f 1 1725.7099 0.0019 
f 1 1883.6007 0.0020 
f 8 1909.9075 0.0025 
f 28 2436.2605 0.0015 
f 35 2620.5189 0.0014 
g 0 1711.9393 0.0071 
g 5 1909.6401 0.0054 
g 18 2423.6690 0.0040 
g 23 2621.3499 0.0027 
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Table 10 

Planetary Parameters of TOI-1136 
Parameter Symbol planet b planet c planet d planet e planet f planet g 
From Transit Modeling 
Planet/Star Radius Ratio | Re 0.0180+1:9е-3 0.02725123е-4 0.04379* 701471 0.024971 660-4 0.03671*52: 4 0.023971:04-3 
Impact Parameter b 0.70573 8-3 0.15514е-! 024514е-! 0.377232 0.421111 0.31711-1 
Scaled Semimajor Axis a/R, 11.202571 14.802 1 2350559 1 30.81+77е-1 38.5679:6-1 50.6* 115100 
Transit Duration (hr) Ti 2.07168 327482 4. ux 445*l3-l 4.961201 5.801201 
Orbital Inclination (deg) i 86.44127е- 89.4243 30-1 89.4173 8¢— 89.312971 во 38222-1 89.65118е—1 
From Stable TTV Solutions! 
Mass Ratio ть/т, 0. 000008853 55-06 0.00001780128— 06 0.00002349*2-2¢— 06 0.00001594*3: x Е 0.00002452111е— 05 0. 0000140475 5:05 
Orbital Period (days) Pow 4.172787 7-36-04 6.25725:17е-04 12.51937*11:-04 18.799211-7е—03 26.3162+17е—03 39.5387+3:6е—03 
Mean Anomaly (deg) M 241} i 681% 120%; UN 170* $2: 100 2353 101 rigen 
Orbital Eccentricity e 0.031734е-02 0:117:226-% 0.0161} 32-93 0.0571} 9¢~ 05 0.012470 0.036740 
Argument of Pericenter (deg) ш 451266101 -113*72:100 118226701 661000 цо 9101" вте 
Eccentricity Vector? ecosw 0.016+0029 —0.046*0020 —0.00917 9.308! 0.022+000% —0.0033*00075 0.0021*00076 
Eccentricity Vector esinw 0.0127 0030 —0.11*0023 -0.0076:0:0146 —0.051*00H —0.0012700127 -0.03410021 
Longitude of Ascending Node (deg) Q 0 (fixed) 0 (fixed) 0 (fixed) 0 (fixed) 0 (fixed) 0 (fixed) 
Orbital Inclination (deg) i 90 (fixed) 90 (fixed) 90 (fixed) 90 (fixed) 90 (fixed) 90 (fixed) 
From RM Modeling 
Sky-projected Stellar Obliquity (deg) À SES 
Stellar Obliquity (deg) 17 <28° (95%) 
Rotational Broadening (km s7!) vsini, 67+ 0.6 
Derived Parameters 
Planetary Radius (Ra) R, 1.907015 2.879 006 4.627 0072 2.6390 3.88 2.53101) 
Planetary Mass (Me) Mp 3.011029 6.0713 8.013 5.4110 8.3158 4.8137 


Note. 1. Reported as osculating Keplerian elements at BJD = 2458680, which is close to the first TESS observation. We note that the osculating orbital periods that depend sensitively on if there had just been a close 


encounter of planets. To compute the orbital period ratio between neighboring planets and the deviation from MMR A, we used the average orbital period from N-body integration in Section 7.1. 2. We also report the 


eccentricity vectors directly because the arguments of pericenter often wrap around 27. 
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